Project overview

The goal of the Pyromania project is to test how terrestrial subsides (plant biomass loading or “browning”) and burning influence aquatic productivity, water quality/chemistry, and trophic transfer. We used a manipulative experiment to assess a range of plant material quantities (0-400g per tank) and fire treatment (burned vs unburned material) and the non-linearity of these effects on aquatic systems through 4 time-point samplings. We used 400L aquatic mesocosms and ran the experiment for ~90d in 2021-2022.

Figure 1. Pyromania experimental setup

DATA SETS
This data set is among 3 to be generated for the project and focuses on:

  1. dissolved organic carbon (DOC) concentrations
  2. total dissolved nitrogen (TDN) concentrations
  3. net ecosystem productivity and respiration (NEP, R) using whole mesocosm dissolved oxygen (O2) measurements
  4. burning effects on donor plant material using elemental analysis
  5. greenhouse gas concentrations (carbon dioxide (CO2), methane (CH4)) from mesocosms
  6. trophic transfer using 15N-labeled sage plants to trace N from the plants into plankton.

TIME POINTS

Time-0 (T0) = before the addition of plant materials, plankton were stocked at this stage
Time-1 (T1) = 10 days after the addition of plant materials
Time-2 (T2) = 31 days after …
Time-3 (T3) = 59 days after …
Time-4 (T4) = 89 days after …

General notes on GAM analyses
We fit the generalized additive models (GAMs) via restricted maximum likelihood (REML) to give stable results with the smoothing parameter (sp) to determine the non-linear relationship between response variables and plant-biomass loading (x-axis). We use automatic smoothing with k value generated automatically from the models, which will set the line ‘wiggliness’. Too low and the relationship becomes linear; too high, and the wiggliness goes haywire.

When using the non-linear smoothing, this is the s(x). When the variable is inside the smooth function, this accounts for the nonlinear shape. We do not use additive non-linear smoothing, which is when two smoothers together, as s(x1) + s(x2), instead we use factor-smooth interaction (detailed below). In addition, we use Treatment (and occasionally plankton size fractions, or Type) as predictors outside of the smooth terms s(x1); this allows for linearity. Continuous variables are rarely linear in GAMs, however, setting categorical variables as linear predictors is more common.

Factor-smooth interactions are written as s(x1 by = fac). This sets different smoothers for different variables of “fac”. Usually, the different smoothers are combined with a varying intercept in case the different categories are different in means and slopes, this would be by adding the fac + s(x1 by = fac), where the +fac allows for a different slope. Similarly, in the absence of by = fac, the smoother is considered a global smoother s(x1), fitting a single line to all the data. If a global smoother is combined with a factor term, then this is akin to varying the intercept but keeping the same slope: fac + s(x1).

The EDF - effective degrees of freedom equate with wiggliness, where edf =1 is a straight line, and higher edfs as more wiggly. GAM smoother significance is described as not being able to draw a horizontal line through the data. Finally, it is also advised to check model concurvity, which is the collinearity with models from 0-1.

DOC

Import the data for DOC and TDN and do a loop to clean up all files and make stacked data in single df. This will take the raw data files, align metadata, and filter to make a new df for models.

detach("package:dplyr", unload = TRUE)
library(dplyr)

## import treatment IDs
IDs<-read.csv("data/treatment.IDs.csv")

##### grab files in a list
Total.DOC.files <- list.files(path="data/DOC.TN", pattern = "csv$", full.names = T)

##### what are the file names, sans extensions using package 'tools'
file.names<-file_path_sans_ext(list.files(path="data/DOC.TN", pattern = "csv$", full.names = F))
############ formatting all data in for loop
  for(i in 1:length(Total.DOC.files))
    {
  data<-read.csv(Total.DOC.files[i], sep=",")
  data<-data[,c(1,3,4)] # removed columns we don't need
  data$File<-Total.DOC.files[i]
  colnames(data)<-c("Tank", "DOC..mg.L", "TN..mg.L", "File")
  data$Tank<- IDs$Tank
  data$Tank<-as.numeric(as.character(data$Tank)) # make the column of samples all numeric
  data <- data[!is.na(as.numeric(as.character(data$Tank))),] # remove all rows that aren't numeric/tanks
  data$Treatment<-IDs$Treatment
  data$plant.mass..g<-IDs$plant.mass..g
  make.names(assign(paste(file.names[i], sep=""), data)) # make the file name the name of new df for loop df
  }
########## this is the end of the loop

#see all dfs you've made, the above will be df matching their file names
# ls() 

#Combine files from loop to single df
DOC.df<-rbind(DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4)

DOC.df$File <- sapply(strsplit(DOC.df$File, "/"), `[`, 3) # extract sample names

# alternative way to code the above
#give the 10th-24th character of the file name, removing the rest
#DOC.df$File<-substr(DOC.df$File, 13, 27) 

#alternatively
# remove the 9 letters ('^.) at start 
# remove the 4 letters (.$') at end
#DOC.df$File<-gsub('^.........|....$', '', DOC.df$File) 

# if equals DOC_T0_11052021 then, T0, if not then T1
DOC.df$Time.point<- as.factor(ifelse(DOC.df$File=="DOC_T0.csv", "T0",
         ifelse (DOC.df$File=="DOC_T1.csv", "T1",
           ifelse (DOC.df$File=="DOC_T2.csv", "T2", 
                   ifelse(DOC.df$File=="DOC_T3.csv", "T3", "T4")))))

#rearrange
DOC.df<- DOC.df %>% 
  select(File, Time.point, Treatment, Tank, plant.mass..g, DOC..mg.L, TN..mg.L) 

DOC.df$Treatment<-as.factor(DOC.df$Treatment)

Analyze DOC at each time point. Run model selection and produce plots for each individual timepoint, later pooled into a 5 panel figure.


######## T0 model
m1.DOC.T0<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")

m2.DOC.T0<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

m3.DOC.T0<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

T0.DOC.AIC<-AIC(m1.DOC.T0, m2.DOC.T0, m3.DOC.T0)
# best is smoother solo

summary(m3.DOC.T0)
anova.gam(m3.DOC.T0)
gam.check(m3.DOC.T0, rep=1000)
draw(m3.DOC.T0)
concrvity(m3.DOC.T0)
par(mfrow = c(2, 2))
plot(m3.DOC.T0, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T0<-plot_difference(
  m1.DOC.T0,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T0.mod.plot<-
  plot_smooths(
  model = m3.DOC.T0,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-0") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting


######## T1 model
m1.DOC.T1<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")

m2.DOC.T1<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

m3.DOC.T1<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

T1.DOC.AIC<-AIC(m1.DOC.T1, m2.DOC.T1, m3.DOC.T1)
# best is smooth by factor

summary(m1.DOC.T1)
anova.gam(m1.DOC.T1)
gam.check(m1.DOC.T1, rep=1000)
draw(m1.DOC.T1)
concrvity(m1.DOC.T1)
par(mfrow = c(2, 2))
plot(m1.DOC.T1, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T1<-plot_difference(
  m1.DOC.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T1.mod.plot<-
  plot_smooths(
  model = m1.DOC.T1,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-10") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting


# effect of treatment, smoothing significant across both treatments
# DOC higher in  unburned, relative to burned


########## T2
m1.DOC.T2<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")

m2.DOC.T2<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

m3.DOC.T2<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

T2.DOC.AIC<-AIC(m1.DOC.T2, m2.DOC.T2, m3.DOC.T2)

# best is smooth by factor

summary(m1.DOC.T2)
anova.gam(m1.DOC.T2)
gam.check(m1.DOC.T2, rep=1000)
draw(m1.DOC.T2)
concrvity(m1.DOC.T2)
par(mfrow = c(2, 2))
plot(m1.DOC.T2, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T2<-plot_difference(
  m1.DOC.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T2.mod.plot<-
  plot_smooths(
  model = m1.DOC.T2,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-31") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

# NO effect of treatment, smoothing significant across both treatments
# DOC equivalent in burned and unburned 
# DOC more variable/wonky across gradient in burned


########## T3
m1.DOC.T3<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")

m2.DOC.T3<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m3.DOC.T3<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

T3.DOC.AIC<-AIC(m1.DOC.T3, m2.DOC.T3, m3.DOC.T3)
# best by factor smooth

summary(m1.DOC.T3)
anova.gam(m1.DOC.T3)
gam.check(m1.DOC.T3, rep=1000)
draw(m1.DOC.T3)
concrvity(m1.DOC.T3)
par(mfrow = c(2, 2))
plot(m1.DOC.T3, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T3<-plot_difference(
  m1.DOC.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T3.mod.plot<-
  plot_smooths(
  model = m1.DOC.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-59") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

# effect of treatment, smoothing significant across both treatments
# DOC higher in  burned vs. unburned


########## T4
m1.DOC.T4<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")

m2.DOC.T4<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

m3.DOC.T4<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

T4.DOC.AIC<-AIC(m1.DOC.T4, m2.DOC.T4, m3.DOC.T4)

# best is global
summary(m3.DOC.T4)
anova.gam(m3.DOC.T4)
gam.check(m3.DOC.T4, rep=1000)
draw(m3.DOC.T4)
concrvity(m3.DOC.T4)
par(mfrow = c(2, 2))
plot(m3.DOC.T4, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T4<-plot_difference(
  m1.DOC.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T4.mod.plot<-
  plot_smooths(
  model = m3.DOC.T4,
  series = plant.mass..g
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-89") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

# no effect of treatment, smoothing significant across both treatments
# DOC equivalent in  burned and unburned


mod.rep<-rep(c("~Treatment + s(plant.mass..g, by= Treatment)", 
              "~Treatment + s(plant.mass..g)", 
              "~s(plant.mass..g)"), times=5)

mod.DOC.df<- data.frame(mod.rep)

AIC.DOC<-bind_rows(T0.DOC.AIC, T1.DOC.AIC, T2.DOC.AIC, T3.DOC.AIC, T4.DOC.AIC)
AIC.DOC.mod<-cbind(mod.DOC.df, AIC.DOC)

write.csv(AIC.DOC.mod, "output/AIC models/AIC.DOC.csv")

DOC tables

Table: Results for DOC Time-0

anova.gam(m3.DOC.T0)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                   edf Ref.df     F p-value
## s(plant.mass..g) 1.04   1.08 1.341    0.24

Table: Results for DOC Time-1

anova.gam(m1.DOC.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.035   0.853
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   6.482  7.532 34.39  <2e-16
## s(plant.mass..g):Treatmentunburned 1.568  1.929 59.34  <2e-16

Table: Results for DOC Time-2

anova.gam(m1.DOC.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.035   0.853
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   6.482  7.532 34.39  <2e-16
## s(plant.mass..g):Treatmentunburned 1.568  1.929 59.34  <2e-16

Table: Results for DOC Time-3

anova.gam(m1.DOC.T3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 12.32 0.00182
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.051  2.532 94.00  <2e-16
## s(plant.mass..g):Treatmentunburned 2.202  2.714 56.55  <2e-16

Table: Results for DOC Time-4

anova.gam(m3.DOC.T4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F p-value
## s(plant.mass..g) 1.928  2.385 29.8  <2e-16

DOC figures

Compile raw plots and model-diff plots for final figures.

###### compile the plots with effect plots

extract.legend <- get_legend(
  # create some space to the left of the legend
  DOC.T1.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))


DOC.alltimes<-plot_grid(
  DOC.T0.mod.plot+ theme(legend.position = "none"),
  DOC.T1.mod.plot+ theme(legend.position = "none"),
  DOC.T2.mod.plot+ theme(legend.position = "none"),
  DOC.T3.mod.plot+ theme(legend.position = "none"),
  DOC.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,8,3), ncol=6)
ggsave("figures/DOC.alltimes.mod.pdf", height=4, width=15)

DOC.alltimes

Model differences between the two factor-smoothers. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

DOC.mod.diffs<-plot_grid(
  DOC.diff.T0+ theme(legend.position = "none")+ ggtitle("DOC - Day-0"),
  DOC.diff.T1+ theme(legend.position = "none")+ ggtitle("DOC - Day-10"),
  DOC.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  DOC.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  DOC.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8,8), ncol=5)
ggsave("figures/DOC.mod.diffs.pdf", height=4, width=14)

DOC.mod.diffs

TDN

Total N analysis and plots, running models and making model-diff plots.

TN.df<-DOC.df

######## T0 model
m1.TN.T0<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T0", data = TN.df, method = "REML")

m2.TN.T0<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T0", data = TN.df, method = "REML")

m3.TN.T0<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T0", data = TN.df, method = "REML")

T0.TN.AIC<-AIC(m1.TN.T0, m2.TN.T0, m3.TN.T0)
# best smooth by factor

summary(m1.TN.T0)
anova.gam(m1.TN.T0)
gam.check(m1.TN.T0, rep=1000)
draw(m1.TN.T0)
concrvity(m1.TN.T0)
par(mfrow = c(2, 2))
plot(m1.TN.T0, all.terms = TRUE, page=1)

# model predictions
TN.diff.T0<-plot_difference(
  m1.TN.T0,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T0.mod.plot<-
  plot_smooths(
  model = m1.TN.T0,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-0") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting

# no treatment effect, compare to simplified model (p=0.901)

######## T1 model
m1.TN.T1<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T1", data = TN.df, method = "REML")

m2.TN.T1<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T1", data = TN.df, method = "REML")

m3.TN.T1<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T1", data = TN.df, method = "REML")

T1.TN.AIC<-AIC(m1.TN.T1, m2.TN.T1, m3.TN.T1)
#best global only

summary(m3.TN.T1)
anova.gam(m3.TN.T1)
gam.check(m3.TN.T1, rep=1000)
draw(m3.TN.T1)
concrvity(m3.TN.T1)
par(mfrow = c(2, 2))
plot(m3.TN.T1, all.terms = TRUE, page=1)

# model predictions
TN.diff.T1<-plot_difference(
  m1.TN.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T1.mod.plot<-
  plot_smooths(
  model = m3.TN.T1,
  series = plant.mass..g
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-10") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting

# TN smoother significant for unburned but not burned (p=0.007)


######## T2 model
m1.TN.T2<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T2", data = TN.df, method = "REML")

m2.TN.T2<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T2", data = TN.df, method = "REML")

m3.TN.T2<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T2", data = TN.df, method = "REML")

T2.TN.AIC<-AIC(m1.TN.T2, m2.TN.T2, m3.TN.T2)
#best global with treatment term

summary(m2.TN.T2)
anova.gam(m2.TN.T2)
gam.check(m2.TN.T2, rep=1000)
draw(m2.TN.T2)
concrvity(m2.TN.T2)
par(mfrow = c(2, 2))
plot(m2.TN.T2, all.terms = TRUE, page=1)

# model predictions
TN.diff.T2<-plot_difference(
  m1.TN.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T2.mod.plot<-
  plot_smooths(
  model = m2.TN.T2,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-31") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting

# Near treatment effect p=0.053, higher TN in unburned
# smoother signif: 0.042


######## T3 model
m1.TN.T3<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T3", data = TN.df, method = "REML")

m2.TN.T3<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T3", data = TN.df, method = "REML")

m3.TN.T3<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T3", data = TN.df, method = "REML")

T3.TN.AIC<-AIC(m1.TN.T3, m2.TN.T3, m3.TN.T3)
#best with smooth by factor term

summary(m1.TN.T3)
anova.gam(m1.TN.T3)
gam.check(m1.TN.T3, rep=1000)
draw(m1.TN.T3)
concrvity(m1.TN.T3)
par(mfrow = c(2, 2))
plot(m1.TN.T3, all.terms = TRUE, page=1)

# model predictions
TN.diff.T3<-plot_difference(
  m1.TN.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T3.mod.plot<-
  plot_smooths(
  model = m1.TN.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-59") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting

# Near treatment effect p=0.075, trend for higher TN in burned
# smoother signif: at <0.001 for both


######## T4 model
m1.TN.T4<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T4", data = TN.df, method = "REML")

m2.TN.T4<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T4", data = TN.df, method = "REML")

m3.TN.T4<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T4", data = TN.df, method = "REML")

T4.TN.AIC<-AIC(m1.TN.T4, m2.TN.T4, m3.TN.T4)
#best  with smooth by factor term

summary(m1.TN.T4)
anova.gam(m1.TN.T4)
gam.check(m1.TN.T4, rep=1000)
draw(m1.TN.T4)
concrvity(m1.TN.T4)
par(mfrow = c(2, 2))
plot(m1.TN.T4, all.terms = TRUE, page=1)

# model predictions
TN.diff.T4<-plot_difference(
  m1.TN.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T4.mod.plot<-
  plot_smooths(
  model = m1.TN.T4,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-89") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting


# effect of treatment (higher TN in the burned) (p=0.020)
# significant smoother effect for burned treatment only (p=0.032)

mod.rep<-rep(c("~Treatment + s(plant.mass..g, by= Treatment)", 
              "~Treatment + s(plant.mass..g)", 
              "~s(plant.mass..g)"), times=5)

mod.TN.df<- data.frame(mod.rep)

AIC.TN<-bind_rows(T0.TN.AIC, T1.TN.AIC, T2.TN.AIC, T3.TN.AIC, T4.TN.AIC)
AIC.TN.mod<-cbind(mod.TN.df, AIC.TN)
write.csv(AIC.TN.mod, "output/AIC models/AIC.TN.mod.csv")

TDN tables

Results for TDN Time-0

anova.gam(m1.TN.T0)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.879   0.357
## 
## Approximate significance of smooth terms:
##                                    edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned     1      1 0.009  0.9266
## s(plant.mass..g):Treatmentunburned   1      1 6.303  0.0186

Table: Results for TDN Time-1

anova.gam(m3.TN.T1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F p-value
## s(plant.mass..g) 2.848  3.492 5.72 0.00274

Table: Results for TDN Time-2

anova.gam(m2.TN.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ Treatment + s(plant.mass..g)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 4.122  0.0532
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F p-value
## s(plant.mass..g) 3.207  3.921 4.87 0.00425

Table: Results for TDN Time-3

anova.gam(m1.TN.T3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df   F p-value
## Treatment  1 3.5  0.0752
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value
## s(plant.mass..g):Treatmentburned   4.359  5.269 23.03  < 2e-16
## s(plant.mass..g):Treatmentunburned 2.457  3.022 10.52 0.000194

Table: Results for TDN Time-4

anova.gam(m1.TN.T4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 6.231  0.0196
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.417  2.973 3.613  0.0318
## s(plant.mass..g):Treatmentunburned 1.000  1.000 0.531  0.4731

TDN figures

Compile raw plots and model-diff plots for final figures.

###### compile the plots with effect plots

TN.mod.alltimes<-plot_grid(
  TN.T0.mod.plot+ theme(legend.position = "none"), 
  TN.T1.mod.plot+ theme(legend.position = "none"),
  TN.T2.mod.plot+ theme(legend.position = "none"),
  TN.T3.mod.plot+ theme(legend.position = "none"),
  TN.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,8,3), ncol=6)
ggsave("figures/TN.mods.plots.pdf", height=4, width=13)

TN.mod.alltimes

Model differences between the two factor-smoothers. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

TN.mod.diffs<-plot_grid(
  TN.diff.T0+ theme(legend.position = "none")+ ggtitle("Day-0"),
  TN.diff.T1+ theme(legend.position = "none")+ ggtitle("Day-10"),
  TN.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  TN.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  TN.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8,8,3), ncol=6)
ggsave("figures/TN.mod.diffs.pdf", height=4, width=13)

TN.mod.diffs

YSI

Import YSI data and produce plots of changes in O2% and net ecosystem productivity (NEP) and respiration (R). The YSI data includes Temp, pH, dissolved oxygen (percent and concentration), and conductivity. Here, we will pull in the raw data and make the new metris NEP and R, determined from differences in DO% from dawn-dusk (NEP) and dusk-dawn (R) over a 24h period in each time period.

#load YSI data
YSI<-read.csv("data/Pyro_YSI.csv")

# fix date
YSI$Date<-as.character(YSI$Date)
YSI$Date<-as.POSIXct(YSI$Date, format="%m/%d/%Y")
YSI$Date<-as.Date(YSI$Date, format="%m/%d/%Y")


####### Time 1 change in O2 ################

#separate time points
YSI.T1<- YSI[(YSI$Time.point=="T1"),]

#calculate NEP for T1
T1.Prod<-YSI.T1[(YSI.T1$Date == "2021-11-15"),] # dawn and dusk for 12h period
T1.Dawn1<-T1.Prod[(T1.Prod$Dawn..Dusk == "dawn"),] # dawn-1 measurements
T1.Dusk<-T1.Prod[(T1.Prod$Dawn..Dusk == "dusk"),] # dusk measurements

T1.Dawn2<-YSI.T1[(YSI.T1$Date == "2021-11-16"),] # dawn-2 measurements, following AM

# make new dataframe
T1.O2<-(T1.Dawn1[,c(2,4:6)]) 
T1.O2$dawn1<-T1.Dawn1$DO.percent
T1.O2$dusk1<-T1.Dusk$DO.percent
T1.O2$dawn2<-T1.Dawn2$DO.percent

# R = dusk - dawn (PM to AM, O2 change of day 1)
# NEP = dusk - dawn (PM to AM, O2 change of day 2)

T1.O2<- mutate(T1.O2, 
                NEP=dusk1 - dawn1,
                R=dawn2 - dusk1) 

#sort
T1.O2<-T1.O2 %>% 
  arrange(Treatment, plant.mass..g) 


################ ################ ################
####### Time 2 change in O2 ################

#separate time points
YSI.T2<- YSI[(YSI$Time.point=="T2"),]

#calculate NEP for T2
T2.Prod<-YSI.T2[(YSI.T2$Date == "2021-12-06"),] # dawn and dusk for 12h period
T2.Dawn1<-T2.Prod[(T2.Prod$Dawn..Dusk == "dawn"),] # dawn-1 measurements
T2.Dusk<-T2.Prod[(T2.Prod$Dawn..Dusk == "dusk"),] # dusk measurements

T2.Dawn2<-YSI.T2[(YSI.T2$Date == "2021-12-07"),] # dawn-2 measurements, following AM

# make new dataframe
T2.O2<-(T2.Dawn1[,c(2,4:6)]) 
T2.O2$dawn1<-T2.Dawn1$DO.percent
T2.O2$dusk1<-T2.Dusk$DO.percent
T2.O2$dawn2<-T2.Dawn2$DO.percent

T2.O2<- mutate(T2.O2, 
                NEP=dusk1 - dawn1,
                R=dawn2 - dusk1) 

#sort
T2.O2<-T2.O2 %>% 
  arrange(Treatment, plant.mass..g) 


################ ################ ################
####### Time 3 change in O2 ################

#separate time points
YSI.T3<- YSI[(YSI$Time.point=="T3"),]

#calculate NEP for T3
T3.Prod<-YSI.T3[(YSI.T3$Date == "2022-01-03"),] # dawn and dusk for 12h period
T3.Dawn1<-T3.Prod[(T3.Prod$Dawn..Dusk == "dawn"),] # dawn-1 measurements
T3.Dusk<-T3.Prod[(T3.Prod$Dawn..Dusk == "dusk"),] # dusk measurements

T3.Dawn2<-YSI.T3[(YSI.T3$Date == "2022-01-04"),] # dawn-2 measurements, following AM

# make new dataframe
T3.O2<-(T3.Dawn1[,c(2,4:6)]) 
T3.O2$dawn1<-T3.Dawn1$DO.percent
T3.O2$dusk1<-T3.Dusk$DO.percent
T3.O2$dawn2<-T3.Dawn2$DO.percent

T3.O2<- mutate(T3.O2, 
                NEP=dusk1 - dawn1,
                R=dawn2 - dusk1) 

#sort
T3.O2<-T3.O2 %>% 
  arrange(Treatment, plant.mass..g) 

################ ################ ################
####### Time 3 change in O2 ################

#separate time points
YSI.T4<- YSI[(YSI$Time.point=="T4"),]

#calculate NEP for T4
T4.Prod<-YSI.T4[(YSI.T4$Date == "2022-02-02"),] # dawn and dusk for 12h period
T4.Dawn1<-T4.Prod[(T4.Prod$Dawn..Dusk == "dawn"),] # dawn-1 measurements
T4.Dusk<-T4.Prod[(T4.Prod$Dawn..Dusk == "dusk"),] # dusk measurements

T4.Dawn2<-YSI.T4[(YSI.T4$Date == "2022-02-03"),] # dawn-2 measurements, following AM

# make new dataframe
T4.O2<-(T4.Dawn1[,c(2,4:6)]) 
T4.O2$dawn1<-T4.Dawn1$DO.percent
T4.O2$dusk1<-T4.Dusk$DO.percent
T4.O2$dawn2<-T4.Dawn2$DO.percent

T4.O2<- mutate(T4.O2, 
                NEP=dusk1 - dawn1,
                R=dawn2 - dusk1) 

#sort
T4.O2<-T4.O2 %>% 
  arrange(Treatment, plant.mass..g) 

################ ################ ################
# combine T1  T2  T3 T4 timepoints
################ ################ ################
O2.tanks<-rbind(T1.O2,T2.O2, T3.O2, T4.O2)

cols<-c("Time.point", "Treatment", "Tank") # columns to make factors
O2.tanks[cols] <- lapply(O2.tanks[cols], factor) # make all these factors
O2.tanks$plant.mass..g<-as.numeric(O2.tanks$plant.mass..g)

DO and O2%

TIME POINT 1: Change in O2% from dissolved oxygen. First, we will also run model fitting on the raw DO% data to apply the same approach for visualizing changes in oxygen across dawn-dusk-dawn measurements. We will then combine all these plots into multi panel figurs for NEP and R, and DO% for each point of measurement.

#########################################################
##################################################################
# total oxygen % plot for the 3 time points (dawn-dusk-dawn)

m1.dawn1.T1<-gam(dawn1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.dawn1.T1<-gam(dawn1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.dawn1.T1<-gam(dawn1 ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.dawn1.AIC<-AIC(m1.dawn1.T1, m2.dawn1.T1, m3.dawn1.T1)
# global with treatment best


summary(m2.dawn1.T1)
anova.gam(m2.dawn1.T1)
gam.check(m2.dawn1.T1, rep=1000)
draw(m2.dawn1.T1)
concrvity(m2.dawn1.T1)
par(mfrow = c(2, 2))
plot(m2.dawn1.T1, all.terms = TRUE, page=1)
# model predictions
dawn1.diff.T1<-plot_difference(
  m1.dawn1.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dawn1.T1.mod.plot<-
  plot_smooths(
  model = m2.dawn1.T1,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=dawn1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T1.dawn1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting
  
# treatment (p=0.0279) and smoothers significant (p<0.001)

  
####### #### Dusk 1 
m1.dusk1.T1<-gam(dusk1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.dusk1.T1<-gam(dusk1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.dusk1.T1<-gam(dusk1 ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.dusk1.AIC<-AIC(m1.dusk1.T1, m2.dusk1.T1, m3.dusk1.T1)
# model with treatment and global smooth best

summary(m2.dusk1.T1)
anova.gam(m2.dusk1.T1)
gam.check(m2.dusk1.T1, rep=1000)
draw(m2.dusk1.T1)
concrvity(m2.dusk1.T1)
par(mfrow = c(2, 2))
plot(m2.dusk1.T1, all.terms = TRUE, page=1)
# model predictions
dusk1.diff.T1<-plot_difference(
  m1.dusk1.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dusk1.T1.mod.plot<-
  plot_smooths(
  model = m2.dusk1.T1,
  series = plant.mass..g,
  comparison= Treatment
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=dusk1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T1.dusk1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# smoother significant for both treatments 


####### #### Dawn 2
m1.dawn2.T1<-gam(dawn2 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.dawn2.T1<-gam(dawn2 ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.dawn2.T1<-gam(dawn2 ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.dawn2.AIC<-AIC(m1.dawn2.T1, m2.dawn2.T1, m3.dawn2.T1)
# treatment and global smooth best

summary(m2.dawn2.T1)
anova.gam(m2.dawn2.T1)
gam.check(m2.dawn2.T1, rep=1000)
draw(m2.dawn2.T1)
concrvity(m2.dawn2.T1)
par(mfrow = c(2, 2))
plot(m2.dawn2.T1, all.terms = TRUE, page=1)
# model predictions
dawn2.diff.T1<-plot_difference(
  m1.dawn2.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned")),
)

###########  
#plot for the model output on rawdata
dawn2.T1.mod.plot<-
  plot_smooths(
  model = m2.dawn2.T1,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=dawn2, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  
  ggtitle("T1.dawn2")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# smoother significant for both (p<0.001)

#### group plots
O2.T1<-plot_grid(
  dawn1.T1.mod.plot+ theme(legend.position = "none"), 
  dusk1.T1.mod.plot+ theme(legend.position = "none"),
  dawn2.T1.mod.plot+ theme(legend.position = "none"),
  extract.legend, 
  rel_widths = c(8,8,8,3), ncol=4)

TIME POINT 2: Change in O2% from dissolved oxygen

############################################################
##############################################################################
# total oxygen % plot for the 3 time points (dawn-dusk-dawn)

m1.dawn1.T2<-gam(dawn1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.dawn1.T2<-gam(dawn1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.dawn1.T2<-gam(dawn1 ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.dawn1.AIC<-AIC(m1.dawn1.T2, m2.dawn1.T2, m3.dawn1.T2)
# factor by smooth best


summary(m1.dawn1.T2)
anova.gam(m1.dawn1.T2)
gam.check(m1.dawn1.T2, rep=1000)
draw(m1.dawn1.T2)
concrvity(m1.dawn1.T2)
par(mfrow = c(2, 2))
plot(m1.dawn1.T2, all.terms = TRUE, page=1)

# model predictions
dawn1.diff.T2<-plot_difference(
  m1.dawn1.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dawn1.T2.mod.plot<-
  plot_smooths(
  model = m1.dawn1.T2,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=dawn1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T2.dawn1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting
  
# smoothers significant (p<0.001)

  
####### #### Dusk 1 
m1.dusk1.T2<-gam(dusk1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.dusk1.T2<-gam(dusk1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.dusk1.T2<-gam(dusk1 ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.dusk1.AIC<-AIC(m1.dusk1.T2, m2.dusk1.T2, m3.dusk1.T2)
# model with smooth by factor best

summary(m1.dusk1.T2)
anova.gam(m1.dusk1.T2)
gam.check(m1.dusk1.T2, rep=1000)
draw(m1.dusk1.T2)
concrvity(m1.dusk1.T2)
par(mfrow = c(2, 2))
plot(m1.dusk1.T2, all.terms = TRUE, page=1)

# model predictions
dusk1.diff.T2<-plot_difference(
  m1.dusk1.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dusk1.T2.mod.plot<-
  plot_smooths(
  model = m1.dusk1.T2,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=dusk1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T2.dusk1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# smoother significant for both treatments 


####### #### Dawn 2
m1.dawn2.T2<-gam(dawn2 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.dawn2.T2<-gam(dawn2 ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.dawn2.T2<-gam(dawn2 ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.dawn2.AIC<-AIC(m1.dawn2.T2, m2.dawn2.T2, m3.dawn2.T2)
# smooth by factor best

summary(m1.dawn2.T2)
anova.gam(m1.dawn2.T2)
gam.check(m1.dawn2.T2, rep=1000)
draw(m1.dawn2.T2)
concrvity(m1.dawn2.T2)
par(mfrow = c(2, 2))
plot(m1.dawn2.T2, all.terms = TRUE, page=1)

# model predictions
dawn2.diff.T2<-plot_difference(
  m1.dawn2.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned")),
)

###########  
#plot for the model output on rawdata
dawn2.T2.mod.plot<-
  plot_smooths(
  model = m1.dawn2.T2,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=dawn2, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T2.dawn2")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# smoother significant for both (p<0.001)

#### group plots
O2.T2<-plot_grid(
  dawn1.T2.mod.plot+ theme(legend.position = "none"), 
  dusk1.T2.mod.plot+ theme(legend.position = "none"),
  dawn2.T2.mod.plot+ theme(legend.position = "none"),
  extract.legend, 
  rel_widths = c(8,8,8,3), ncol=4)

TIME POINT 3: Change in O2% from dissolved oxygen

############################################################
##############################################################################
#### Dawn1
m1.dawn1.T3<-gam(dawn1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.dawn1.T3<-gam(dawn1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.dawn1.T3<-gam(dawn1 ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.dawn1.AIC<-AIC(m1.dawn1.T3, m2.dawn1.T3, m3.dawn1.T3)
# factor by smooth best


summary(m1.dawn1.T3)
anova.gam(m1.dawn1.T3)
gam.check(m1.dawn1.T3, rep=1000)
draw(m1.dawn1.T3)
concrvity(m1.dawn1.T3)
par(mfrow = c(2, 2))
plot(m1.dawn1.T3, all.terms = TRUE, page=1)

# model predictions
dawn1.diff.T3<-plot_difference(
  m1.dawn1.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dawn1.T3.mod.plot<-
  plot_smooths(
  model = m1.dawn1.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=dawn1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T3.dawn1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting
  
# smoothers significant (p<0.001)

  
####### #### Dusk 1 
m1.dusk1.T3<-gam(dusk1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.dusk1.T3<-gam(dusk1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.dusk1.T3<-gam(dusk1 ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.dusk1.AIC<-AIC(m1.dusk1.T3, m2.dusk1.T3, m3.dusk1.T3)
# model with smooth by factor best

summary(m1.dusk1.T3)
anova.gam(m1.dusk1.T3)
gam.check(m1.dusk1.T3, rep=1000)
draw(m1.dusk1.T3)
concrvity(m1.dusk1.T3)
par(mfrow = c(2, 2))
plot(m1.dusk1.T3, all.terms = TRUE, page=1)

# model predictions
dusk1.diff.T3<-plot_difference(
  m1.dusk1.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dusk1.T3.mod.plot<-
  plot_smooths(
  model = m1.dusk1.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=dusk1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T3.dusk1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect
# smoother significant for burned 


####### #### Dawn 2
m1.dawn2.T3<-gam(dawn2 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.dawn2.T3<-gam(dawn2 ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.dawn2.T3<-gam(dawn2 ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.dawn2.AIC<-AIC(m1.dawn2.T3, m2.dawn2.T3, m3.dawn2.T3)
# smooth by factor best

summary(m1.dawn2.T3)
anova.gam(m1.dawn2.T3)
gam.check(m1.dawn2.T3, rep=1000)
draw(m1.dawn2.T3)
concrvity(m1.dawn2.T3)
par(mfrow = c(2, 2))
plot(m1.dawn2.T3, all.terms = TRUE, page=1)

# model predictions
dawn2.diff.T3<-plot_difference(
  m1.dawn2.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned")),
)

###########  
#plot for the model output on rawdata
dawn2.T3.mod.plot<-
  plot_smooths(
  model = m1.dawn2.T3,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=dawn2, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T3.dawn2")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# global smoother significant(p<0.001)

#### group plots
O2.T3<-plot_grid(
  dawn1.T3.mod.plot+ theme(legend.position = "none"), 
  dusk1.T3.mod.plot+ theme(legend.position = "none"),
  dawn2.T3.mod.plot+ theme(legend.position = "none"),
  extract.legend, 
  rel_widths = c(8,8,8,3), ncol=4)

TIME POINT 4: Change in O2% from dissolved oxygen

############################################################
##############################################################################
# total oxygen % plot for the 3 time points (dawn-dusk-dawn)

#### Dawn1
m1.dawn1.T4<-gam(dawn1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.dawn1.T4<-gam(dawn1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.dawn1.T4<-gam(dawn1 ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.dawn1.AIC<-AIC(m1.dawn1.T4, m2.dawn1.T4, m3.dawn1.T4)
# model with global best


summary(m3.dawn1.T4)
anova.gam(m3.dawn1.T4)
gam.check(m3.dawn1.T4, rep=1000)
draw(m3.dawn1.T4)
concrvity(m3.dawn1.T4)
par(mfrow = c(2, 2))
plot(m3.dawn1.T4, all.terms = TRUE, page=1)

# model predictions
dawn1.diff.T4<-plot_difference(
  m1.dawn1.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dawn1.T4.mod.plot<-
  plot_smooths(
  model = m3.dawn1.T4,
  series = plant.mass..g,
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=dawn1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T4.dawn1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting
  
# smoothers significant (p<0.001)

  
####### #### Dusk 1 
m1.dusk1.T4<-gam(dusk1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.dusk1.T4<-gam(dusk1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.dusk1.T4<-gam(dusk1 ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.dusk1.AIC<-AIC(m1.dusk1.T4, m2.dusk1.T4, m3.dusk1.T4)
# model with smooth by factor best

summary(m1.dusk1.T4)
anova.gam(m1.dusk1.T4)
gam.check(m1.dusk1.T4, rep=1000)
draw(m1.dusk1.T4)
concrvity(m1.dusk1.T4)
par(mfrow = c(2, 2))
plot(m1.dusk1.T4, all.terms = TRUE, page=1)

# model predictions
dusk1.diff.T4<-plot_difference(
  m1.dusk1.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dusk1.T4.mod.plot<-
  plot_smooths(
  model = m1.dusk1.T4,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=dusk1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T4.dusk1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect
# smoother significant for burned 


####### #### Dawn 2
m1.dawn2.T4<-gam(dawn2 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.dawn2.T4<-gam(dawn2 ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.dawn2.T4<-gam(dawn2 ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.dawn2.AIC<-AIC(m1.dawn2.T4, m2.dawn2.T4, m3.dawn2.T4)
# global smooth best

summary(m3.dawn2.T4)
anova.gam(m3.dawn2.T4)
gam.check(m3.dawn2.T4, rep=1000)
draw(m3.dawn2.T4)
concrvity(m3.dawn2.T4)
par(mfrow = c(2, 2))
plot(m3.dawn2.T4, all.terms = TRUE, page=1)

# model predictions
dawn2.diff.T4<-plot_difference(
  m1.dawn2.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned")),
)

###########  
#plot for the model output on rawdata
dawn2.T4.mod.plot<-
  plot_smooths(
  model = m3.dawn2.T4,
  series = plant.mass..g,
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=dawn2, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T4.dawn2")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# global smoother significant(p<0.001)

#### group plots
O2.T4<-plot_grid(
  dawn1.T4.mod.plot+ theme(legend.position = "none"), 
  dusk1.T4.mod.plot+ theme(legend.position = "none"),
  dawn2.T4.mod.plot+ theme(legend.position = "none"),
  extract.legend, 
  rel_widths = c(8,8,8,3), ncol=4)

Combine and export all the O2 data with plot-difference and model AIC tables

#### model differences
O2.mod.diffs<-plot_grid(
  dawn1.diff.T1+ theme(legend.position = "none")+ ggtitle("T1-Dawn1"),
  dusk1.diff.T1+ theme(legend.position = "none")+ ggtitle("Dusk1"),
  dawn2.diff.T1+ theme(legend.position = "none")+ ggtitle("Dawn2"),
  
  dawn1.diff.T2+ theme(legend.position = "none")+ ggtitle("T2-Dawn1"),
  dusk1.diff.T2+ theme(legend.position = "none")+ ggtitle("Dusk1"),
  dawn2.diff.T2+ theme(legend.position = "none")+ ggtitle("Dawn2"),
  
  dawn1.diff.T3+ theme(legend.position = "none")+ ggtitle("T3-Dawn1"),
  dusk1.diff.T3+ theme(legend.position = "none")+ ggtitle("Dusk1"),
  dawn2.diff.T3+ theme(legend.position = "none")+ ggtitle("Dawn2"),
  
  dawn1.diff.T4+ theme(legend.position = "none")+ ggtitle("T4-Dawn1"),
  dusk1.diff.T4+ theme(legend.position = "none")+ ggtitle("Dusk1"),
  dawn2.diff.T4+ theme(legend.position = "none")+ ggtitle("Dawn2"),
  rel_widths = c(8,8,8), ncol=3, nrow=4)

ggsave("figures/O2.mod.diffs.pdf", height=10, width=7)

#### model and raw data
O2.mods<-plot_grid(
  O2.T1+ theme(legend.position = "none")+ ggtitle("Day-10"),
  O2.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  O2.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  O2.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8), ncol=1, nrow=4)

ggsave("figures/O2.mod.pdf", height=11, width=9)


#bind the AIC tables
AIC.O2<-bind_rows(T1.dawn1.AIC, T1.dusk1.AIC, T1.dawn2.AIC,
                  T2.dawn1.AIC, T2.dusk1.AIC, T2.dawn2.AIC,
                  T3.dawn1.AIC, T3.dusk1.AIC, T3.dawn2.AIC,
                  T4.dawn1.AIC, T4.dusk1.AIC, T4.dawn2.AIC)

# make a model column
mod.rep12<-rep(c("~Treatment + s(plant.mass..g, by= Treatment)", 
              "~Treatment + s(plant.mass..g)", 
              "~s(plant.mass..g)"), times=12)

mod.O2.df<- data.frame(mod.rep12)
#bind table

AIC.O2.mod<-cbind(mod.O2.df, AIC.O2)

write.csv(AIC.O2.mod, "output/AIC models/AIC.O2.mod.csv")

NEP and R

Generate dataframes for NEP and R change in O2. We will use NEP and R, run gam model fits, and produce individual figures for each time point.

First, we will use NEP for productivity measurements.

####### Time 1

m1.NEP.T1<-gam(NEP ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.NEP.T1<-gam(NEP ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.NEP.T1<-gam(NEP ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.NEP.AIC<-AIC(m1.NEP.T1, m2.NEP.T1, m3.NEP.T1)
# model with plot smooth by factor not different from reduced model, go with smooth by factor

summary(m2.NEP.T1)
anova.gam(m2.NEP.T1)
gam.check(m2.NEP.T1, rep=1000)
draw(m2.NEP.T1)
concrvity(m2.NEP.T1)
par(mfrow = c(2, 2))
plot(m2.NEP.T1, all.terms = TRUE, page=1)


#### see this https://cran.r-project.org/web/packages/tidymv/vignettes/plot-smooths.html
# The difference smooth is difference between the smooths of two conditions (two levels in a factor). 
# Portions of the difference smooth confidence interval that do not include 0 are shaded in red.

# model predictions
NEP.diff.T1<-plot_difference(
  m2.NEP.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
NEP.T1.mod.plot<-
  plot_smooths(
  model = m2.NEP.T1,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=NEP, color=Treatment)) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-20, 50)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Production (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect (p=0.110), smoothers significant (p<0.006)

####### Time 2

m1.NEP.T2<-gam(NEP ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.NEP.T2<-gam(NEP ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.NEP.T2<-gam(NEP ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.NEP.AIC<-AIC(m1.NEP.T2, m2.NEP.T2, m3.NEP.T2)
# model with plot smooth by factor not different from reduced model, go with smooth by factor


summary(m2.NEP.T2)
anova.gam(m2.NEP.T2)
gam.check(m2.NEP.T2, rep=1000)
draw(m2.NEP.T2)
concrvity(m2.NEP.T2)
par(mfrow = c(2, 2))
plot(m2.NEP.T2, all.terms = TRUE, page=1)

## plot for the model output on rawdata
# model predictions
NEP.diff.T2<-plot_difference(
  m2.NEP.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

NEP.T2.mod.plot<-
  plot_smooths(
  model = m2.NEP.T2,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=NEP, color=Treatment)) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-20, 50)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Production (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# treatment effect (p=0.007)
# smoother significant for burned (p=0.002) but not unburned (p=0.326)


####### Time 3

m1.NEP.T3<-gam(NEP ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.NEP.T3<-gam(NEP ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.NEP.T3<-gam(NEP ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.NEP.AIC<-AIC(m1.NEP.T3, m2.NEP.T3, m3.NEP.T3)
# model with plot smooth by factor not different from reduced model, go with smooth by factor


summary(m2.NEP.T3)
anova.gam(m2.NEP.T3)
gam.check(m2.NEP.T3, rep=1000)
draw(m2.NEP.T3)
concrvity(m2.NEP.T3)
par(mfrow = c(2, 2))
plot(m2.NEP.T3, all.terms = TRUE, page=1)

# model predictions
NEP.diff.T3<-plot_difference(
  m2.NEP.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

#plot for the model output on rawdata
NEP.T3.mod.plot<-
  plot_smooths(
  model = m2.NEP.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=NEP, color=Treatment)) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-20, 50)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Production (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# treatment effect (p=0.009)
# smoother significant for burned (p<0.001) but not unburned (p=0.053)


####### Time 4

m1.NEP.T4<-gam(NEP ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.NEP.T4<-gam(NEP ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.NEP.T4<-gam(NEP ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.NEP.AIC<-AIC(m1.NEP.T4, m2.NEP.T4, m3.NEP.T4)
# model with plot smooth by factor not different from reduced model, go with smooth by factor


summary(m1.NEP.T4)
anova.gam(m1.NEP.T4)
gam.check(m1.NEP.T4, rep=1000)
draw(m1.NEP.T4)
concrvity(m1.NEP.T4)
par(mfrow = c(2, 2))
plot(m1.NEP.T4, all.terms = TRUE, page=1)

# model predictions
NEP.diff.T4<-plot_difference(
  m1.NEP.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

#plot for the model output on rawdata
NEP.T4.mod.plot<-
  plot_smooths(
  model = m1.NEP.T4,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=NEP, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-20, 50)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Production (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect (p=0.118)
# smoother significant for burned (p=0.020) but not unburned (p=0.327)

mod.RNEP<-rep(c("~Treatment + s(plant.mass..g, by= Treatment)", 
              "~Treatment + s(plant.mass..g)", 
              "~s(plant.mass..g)"), times=4)

mod.RNEP.df<- data.frame(mod.RNEP)

AIC.NEP<-bind_rows(T1.NEP.AIC, T2.NEP.AIC, T3.NEP.AIC, T4.NEP.AIC)
AIC.NEP.mod<-cbind(mod.RNEP.df, AIC.NEP)
write.csv(AIC.NEP.mod, "output/AIC models/AIC.NEP.mod.csv")

Table: Results for Time-1 NEP

anova.gam(m1.NEP.T1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## NEP ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 2.756    0.11
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.478  3.046 13.44 2.3e-05
## s(plant.mass..g):Treatmentunburned 1.690  2.090  6.36 0.00585

Table: Results for Time-2 NEP

anova.gam(m1.NEP.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## NEP ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 8.479 0.00745
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   1.967  2.431 7.392 0.00227
## s(plant.mass..g):Treatmentunburned 1.000  1.000 1.002 0.32638

Table: Results for Time-3 NEP

anova.gam(m1.NEP.T3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## NEP ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 8.118 0.00948
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value
## s(plant.mass..g):Treatmentburned   3.414  4.166 8.236 0.000334
## s(plant.mass..g):Treatmentunburned 3.124  3.821 2.940 0.053008

Table: Results for Time-4 NEP

anova.gam(m1.NEP.T4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## NEP ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df    F p-value
## Treatment  1 2.62   0.118
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.757  3.382 3.717  0.0203
## s(plant.mass..g):Treatmentunburned 1.000  1.000 1.002  0.3268

Now, we will go through Respiration models and individual plots.

####### Time 1
m1.R.T1<-gam(R ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.R.T1<-gam(R ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.R.T1<-gam(R ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.R.AIC<-AIC(m1.R.T1, m2.R.T1, m3.R.T1)
# model with global best


summary(m3.R.T1)
anova.gam(m3.R.T1)
gam.check(m3.R.T1, rep=1000)
draw(m3.R.T1)
concrvity(m3.R.T1)
par(mfrow = c(2, 2))
plot(m3.R.T1, all.terms = TRUE, page=1)
# model predictions
R.diff.T1<-plot_difference(
  m1.R.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
R.T1.mod.plot<-
  plot_smooths(
  model = m3.R.T1,
  series = plant.mass..g,
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=R, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-40, 10)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Respiration (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect (p=0.229), smoothers significant (p<0.001)

####### Time 2

m1.R.T2<-gam(R ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.R.T2<-gam(R ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.R.T2<-gam(R ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.R.AIC<-AIC(m1.R.T2, m2.R.T2, m3.R.T2)
# model with global + treatment best

summary(m2.R.T2)
anova.gam(m2.R.T2)
gam.check(m2.R.T2, rep=1000)
draw(m2.R.T2)
concrvity(m2.R.T2)
par(mfrow = c(2, 2))
plot(m2.R.T2, all.terms = TRUE, page=1)
# model predictions
R.diff.T2<-plot_difference(
  m2.R.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
R.T2.mod.plot<-
  plot_smooths(
  model = m2.R.T2,
  series = plant.mass..g, 
  comparison= Treatment,
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=R, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  coord_cartesian(ylim=c(-40, 10)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Respiration (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# slight treatment effect (p=0.085)
# smoother significant for both burned and unburned (p<0.008)


####### Time 3
m1.R.T3<-gam(R ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.R.T3<-gam(R ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.R.T3<-gam(R ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.R.AIC<-AIC(m1.R.T3, m2.R.T3, m3.R.T3)
# model smmoth by factor best


summary(m1.R.T3)
anova.gam(m1.R.T3)
gam.check(m1.R.T3, rep=1000)
draw(m1.R.T3)
concrvity(m1.R.T3)
par(mfrow = c(2, 2))
plot(m1.R.T3, all.terms = TRUE, page=1)
# model predictions
R.diff.T3<-plot_difference(
  m1.R.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
R.T3.mod.plot<-
  plot_smooths(
  model = m1.R.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=R, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-40, 10)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Respiration (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# treatment effect (p=0.027)
# smoother significant for burned and unnburned (p<0.001)


####### Time 4

m1.R.T4<-gam(R ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.R.T4<-gam(R ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.R.T4<-gam(R ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.R.AIC<-AIC(m1.R.T4, m2.R.T4, m3.R.T4)
# model with global + treatment best


summary(m1.R.T4)
anova.gam(m1.R.T4)
gam.check(m1.R.T4, rep=1000)
draw(m1.R.T4)
concrvity(m1.R.T4)
par(mfrow = c(2, 2))
plot(m1.R.T4, all.terms = TRUE, page=1)
# model predictions
R.diff.T4<-plot_difference(
  m1.R.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
R.T4.mod.plot<-
  plot_smooths(
  model = m1.R.T4,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=R, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-40, 10)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Respiration (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect (p=0.078)
# smoother significant for burned (p=0.004) but not unburned (p=0.965)


R.mod.plot<-plot_grid(
  R.T1.mod.plot+ theme(legend.position = "none")+ ggtitle("Day-10"),
  R.T2.mod.plot+ theme(legend.position = "none")+ ggtitle("Day-31"),
  R.T3.mod.plot+ theme(legend.position = "none")+ ggtitle("Day-59"),
  R.T4.mod.plot+ theme(legend.position = "none")+ ggtitle("Day-89"), extract.legend,
  rel_widths = c(8,8,8,8,3), ncol=5)

#ggsave("figures/R.mod.plot.long.pdf", height=6, width=12)

#### model differences
R.mod.diffs<-plot_grid(
  R.diff.T1+ theme(legend.position = "none")+ ggtitle("R-Day-10"),
  R.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  R.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  R.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8), ncol=4)

#ggsave("figures/R.mod.diffs.pdf", height=3, width=10)

AIC.R<-bind_rows(T1.R.AIC, T2.R.AIC, T3.R.AIC, T4.R.AIC)
AIC.R.mod<-cbind(mod.RNEP.df, AIC.R)

write.csv(AIC.R.mod, "output/AIC models/AIC.R.mod.csv")

Table: Results for Time-1 R

anova.gam(m3.R.T1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## R ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F  p-value
## s(plant.mass..g) 2.520  3.097 12.8 2.32e-05

Table: Results for Time-2 R

anova.gam(m2.R.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## R ~ Treatment + s(plant.mass..g)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 6.443  0.0186
## 
## Approximate significance of smooth terms:
##                    edf Ref.df  F  p-value
## s(plant.mass..g) 5.710  6.758 10 1.52e-05

Table: Results for Time-3 R

anova.gam(m1.R.T3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## R ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 5.669  0.0268
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F  p-value
## s(plant.mass..g):Treatmentburned   3.762  4.576 13.144 1.02e-05
## s(plant.mass..g):Treatmentunburned 3.274  4.000  7.775 0.000523

Table: Results for Time-4 R

anova.gam(m1.R.T4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## R ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df    F p-value
## Treatment  1 3.38  0.0784
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.927  3.587 5.293 0.00408
## s(plant.mass..g):Treatmentunburned 1.000  1.000 0.002 0.96531

Compiled NEP-R plots with model fits

NEP.R.alltimes.long<-plot_grid(
  NEP.T1.mod.plot+ theme(legend.position = "none"),
  NEP.T2.mod.plot+ theme(legend.position = "none"),
  NEP.T3.mod.plot+ theme(legend.position = "none"),
  NEP.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  R.T1.mod.plot+ theme(legend.position = "none"),
  R.T2.mod.plot+ theme(legend.position = "none"),
  R.T3.mod.plot+ theme(legend.position = "none"),
  R.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,3, 8,8,8,8,3), ncol=5)
ggsave("figures/NEP.R.alltimes.long.pdf", height=7, width=12)

NEP.R.alltimes.long

Model differences between the two factor-smoothers

NEP.R.mod.diffs<-plot_grid(
  NEP.diff.T1+ theme(legend.position = "none"),
  NEP.diff.T2+ theme(legend.position = "none"),
  NEP.diff.T3+ theme(legend.position = "none"),
  NEP.diff.T4+ theme(legend.position = "none"),
  R.diff.T1+ theme(legend.position = "none"),
  R.diff.T2+ theme(legend.position = "none"),
  R.diff.T3+ theme(legend.position = "none"),
  R.diff.T4+ theme(legend.position = "none"),
  rel_widths = c(8,8,8,8,8,8,8,8), ncol=4)
ggsave("figures/NEP.R.mod.diffs.alt.pdf", height=7, width=12)

NEP.R.mod.diffs

Isotopes

Isotopes and C:N for starting materials used in the experiment and plankton fractions sampled at time 1 and time 2.

#########  Time 1 and Time 2
topes<-read.csv("data/Isotopes/Pyro_isotopes.csv")
topes$C.N <-(topes$Total.C..ug/12)/(topes$Total.N..ug/14) # C mol : N mol

cols<-c("Time.point", "Treatment", "Type", "Tank") # columns to make factors
topes[cols] <- lapply(topes[cols], factor) # make all these factors

##### make data frames
# treatment data df
topes.trt<-topes[(topes$Treatment=="burned" | topes$Treatment=="unburned"),]
topes.trt<-droplevels(topes.trt)
topes.trt$Type<-factor(topes.trt$Type, 
                         levels=c("plankton", "POM"))

# control and start plant materials df
topes.controls<-topes[!(topes$Treatment=="burned" | topes$Treatment=="unburned"),]

Total C vs total N for plankton at Time1 and Time2. The model here is a linear regression

############# C vs N plot

#plot for the model output on rawdata
CvN.mod1.T1<-lm(Total.N..ug~Total.C..ug*Treatment, data=topes.trt[(topes.trt$Time.point=="T1"),])

CvN.mod2.T1<-lm(Total.N..ug~Total.C..ug + Treatment, data=topes.trt[(topes.trt$Time.point=="T1"),])

anova(CvN.mod1.T1,CvN.mod2.T1) # no diff, go with simple model 'CvN.mod2.T1'
## Analysis of Variance Table
## 
## Model 1: Total.N..ug ~ Total.C..ug * Treatment
## Model 2: Total.N..ug ~ Total.C..ug + Treatment
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     56 22978                           
## 2     57 23284 -1   -305.75 0.7451 0.3917
Anova(CvN.mod1.T1, type=3)
## Anova Table (Type III tests)
## 
## Response: Total.N..ug
##                       Sum Sq Df  F value  Pr(>F)    
## (Intercept)             1590  1   3.8744 0.05399 .  
## Total.C..ug           101937  1 248.4283 < 2e-16 ***
## Treatment                 15  1   0.0368 0.84859    
## Total.C..ug:Treatment    306  1   0.7451 0.39170    
## Residuals              22978 56                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
CvN.T1.mod.plot<-
  ggplot(data=topes.trt[(topes.trt$Time.point=="T1"),], 
             aes(x=Total.C..ug, y=Total.N..ug, 
                 color=Treatment, # for the points
                 fill=Treatment, # for the SE
                 linetype=Treatment)) + #for the regression lines
  geom_point(aes(shape=Type))+
  geom_smooth(method = "lm", alpha =0.2, size=0.5)+ 
  scale_color_manual(values = c("brown1", "mediumseagreen"))+
  ggtitle("Time-1") + 
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("< 63"~mu,"m")), 
                                expression(paste("> 63"~mu,"m")))) +
  ylim(c(0,600)) +
  xlim(c(0,3200)) +
  ylab("Total N (ug)") +
  xlab("Total C (ug)") +
  Fig.formatting 


#############################
### Time 2 C vs. N
#plot for the model output on rawdata
CvN.mod1.T2<-lm(Total.N..ug~Total.C..ug*Treatment, data=topes.trt[(topes.trt$Time.point=="T2"),])

CvN.mod2.T2<-lm(Total.N..ug~Total.C..ug + Treatment, data=topes.trt[(topes.trt$Time.point=="T2"),])

anova(CvN.mod1.T2,CvN.mod2.T2) # no diff, go with simple model 'CvN.mod2.T2'
## Analysis of Variance Table
## 
## Model 1: Total.N..ug ~ Total.C..ug * Treatment
## Model 2: Total.N..ug ~ Total.C..ug + Treatment
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     56 14419                           
## 2     57 14732 -1   -312.27 1.2128 0.2755
Anova(CvN.mod2.T2, type=2)
## Anova Table (Type II tests)
## 
## Response: Total.N..ug
##             Sum Sq Df  F value Pr(>F)    
## Total.C..ug 204416  1 790.9323 <2e-16 ***
## Treatment      194  1   0.7492 0.3904    
## Residuals    14732 57                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
CvN.T2.mod.plot<-
  ggplot(data=topes.trt[(topes.trt$Time.point=="T2"),], 
             aes(x=Total.C..ug, y=Total.N..ug, 
                 color=Treatment, # for the points
                 fill=Treatment, # for the SE
                 linetype=Treatment)) + #for the regression lines
  geom_point(aes(shape=Type))+
  geom_smooth(method = "lm", alpha =0.2, size=0.5)+
  scale_color_manual(values = c("brown1", "mediumseagreen"))+
  ggtitle("Time-2") + 
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  ylim(c(0,250)) +
  xlim(c(0,1250)) +
  ylab("Total N (ug)") +
  xlab("Total C (ug)") +
  Fig.formatting 

C:N vs. plant material models and plot, fit with GAMs Test if the C:N was consistent across time, or what factors may influence C:N. This is important to verify assumptions. We are using 15N transer between trophic levels, assuming our measurements of efficiency in nitrogen transfer reflect carbon and eRgy transfer between trophic levels.

############# all plankton T1 nd T2
m1.CN <- gam(C.N ~ Treatment + Type + Time.point +
            s(plant.mass..g, by=Treatment),
            data=topes.trt, method="REML", family="gaussian")

m2.CN <- gam(C.N ~ Treatment + Type + Time.point +
            s(plant.mass..g), 
            data=topes.trt, method="REML", family="gaussian")

m3.CN <- gam(C.N ~
            s(plant.mass..g), 
            data=topes.trt, method="REML", family="gaussian")

m4.CN <- gam(C.N ~ Treatment + Type +
            s(plant.mass..g, by=Time.point),
            data=topes.trt, method="REML", family="gaussian")

m5.CN <- gam(C.N ~ Treatment + Type +
            s(plant.mass..g, by=Treatment),
            data=topes.trt, method="REML", family="gaussian")

#best model here
m6.CN <- gam(C.N ~ Type +
            s(plant.mass..g, by=Treatment),
            data=topes.trt, method="REML", family="gaussian")

m7.CN <- gam(C.N ~ Treatment +
            s(plant.mass..g, by=Type),
            data=topes.trt, method="REML", family="gaussian")

m8.CN <- gam(C.N ~ Type +
            s(plant.mass..g, by=Type),
            data=topes.trt, method="REML", family="gaussian")


AIC.CN<-AIC(m1.CN, m2.CN, m3.CN, m4.CN, m5.CN, m6.CN, m7.CN, m8.CN)
## additive model best fit, but no treatment or type effect

summary(m6.CN)
anova.gam(m6.CN)
gam.check(m6.CN, rep=1000)
draw(m6.CN)
concrvity(m6.CN)
par(mfrow = c(1, 2))
plot(m6.CN, all.terms = TRUE, page=1)
# model for smoothing
msmooth.CN<-gam(C.N ~ Type +
            s(plant.mass..g, by=Type),
            data=topes.trt, method="REML", family="gaussian")

# model predictions
CN.diff<-plot_difference(
  m1.CN,
  series = plant.mass..g,
  difference = list(Type = c("plankton", "POM"))
)


###########  

#plot for the model output on rawdata
CN.mod.plot.timepooled<-
  plot_smooths(
  model = msmooth.CN,
  series = plant.mass..g,
  comparison = Type)  + 
  theme(legend.position = "none") +
  geom_point(data=topes.trt, 
             aes(x=plant.mass..g, y=C.N, color=Type, fill=Type)) +
  scale_fill_manual(values = c("deepskyblue4", "darkseagreen"), guide='none') +
  scale_color_manual(name="Plankton", values = c("deepskyblue4", "darkseagreen"), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  theme(legend.position = "right") +
  ggtitle("Days 10 and 31") + 
  coord_cartesian(ylim=c(0, 20)) +
  ylab("C:N") +
  xlab("plant material (g)") +
  Fig.formatting 

CN.mod.plot.timepooled
ggsave("figures/CN.mod.plot.timepooled.pdf", height=4, width=5, encod="MacRoman")

#####

# linear model approach
CN.all.mod<-lm(C.N~ Treatment+Type+Time.point, data=topes.trt, na.action=na.exclude) 
print(Anova(CN.all.mod, type=2), digits=5)
posthoc<-emmeans(CN.all.mod, ~Type)
multcomp::cld(posthoc, Letters=letters)


# All time C.N boxplot
CNbox.all.time<-ggplot(topes.trt, aes(x=Treatment, y=C.N, fill=Type)) + 
  geom_boxplot() +
  geom_point(pch = 21, position = position_jitterdodge(), alpha=0.6) +
  scale_fill_manual(name="Plankton", values = c("deepskyblue4", "darkseagreen"), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  coord_cartesian(ylim=c(0, 20)) +
  ylab("C:N") +
  Fig.formatting

CNbox.all.time
ggsave("figures/CNbox.all.time.pdf", height=4, width=5, encod="MacRoman")
############# all plankton T1
m1.T1.CN <- gam(C.N ~ Treatment + Type +
            s(plant.mass..g, by=Treatment),
        subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m2.T1.CN <- gam(C.N ~ Treatment + Type +
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m3.T1.CN <- gam(C.N ~
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")


AIC.CN.T1<-AIC(m1.T1.CN, m2.T1.CN, m3.T1.CN)
## additive model best fit, but no treatment or type effect

summary(m1.T1.CN)
anova.gam(m1.T1.CN)
gam.check(m1.T1.CN, rep=1000)
draw(m1.T1.CN)
concrvity(m1.T1.CN)
par(mfrow = c(1, 2))
plot(m1.T1.CN, all.terms = TRUE, page=1)
# model for smoothing
msmooth.T1.CN<- gam(C.N ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

# model predictions
CN.diff.T1<-plot_difference(
  m1.T1.CN,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)


###########  

#plot for the model output on rawdata
CN.T1.mod.plot<-
  plot_smooths(
  model = m1.T1.CN,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=C.N, color=Treatment, shape=Type)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  theme(legend.position = "right") +
  ggtitle("Time-1") + 
  coord_cartesian(ylim=c(0, 20)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  ylab("C:N") +
  xlab("plant material (g)") +
  Fig.formatting 

# no effect of type or treatment
# smoother significant for both

############# all plankton T2
m1.T2.CN <- gam(C.N ~ Treatment + Type +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m2.T2.CN <- gam(C.N ~ Treatment + Type +
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m3.T2.CN <- gam(C.N ~
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")


AIC.CN.T2<-AIC(m1.T2.CN, m2.T2.CN, m3.T2.CN)
## additive model best fit, but no treatment or type effect

summary(m1.T2.CN)
anova.gam(m1.T2.CN)
gam.check(m2.T2.CN, rep=1000)
draw(m1.T2.CN)
concrvity(m1.T2.CN)
par(mfrow = c(1, 2))
plot(m1.T2.CN, all.terms = TRUE, page=1)
 ###########  
# model for smoothing
msmooth.T2.CN<- gam(C.N ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

# model predictions
CN.diff.T2<-plot_difference(
  msmooth.T2.CN,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

#plot for the model output on rawdata
CN.T2.mod.plot<-
  plot_smooths(
  model = msmooth.T2.CN,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=C.N, color=Treatment, shape=Type)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  theme(legend.position = "right") +
  ggtitle("Time-2") + 
  coord_cartesian(ylim=c(0, 20)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  ylab("C:N") +
  xlab("plant material (g)") +
  Fig.formatting 

# no effect of treatment, but type important (higher in POM)
# smoother significant for both

Compile CvsN and C:N vs. plant material plots and models

## AIC table
mod.CN<-rep(c("Treatment + Type + s(plant.mass..g, by=Treatment)", 
              "Treatment + Type + s(plant.mass..g)",
              "s(plant.mass..g)"), times=2)

mod.CN.df<- data.frame(mod.CN)
AIC.CN<-bind_rows(AIC.CN.T1, AIC.CN.T2)

AIC.CN.mod<-cbind(mod.CN.df, AIC.CN)
write.csv(AIC.CN.mod, "output/AIC models/AIC.CN.mod.csv")


### get legend
extract.legend.C.N <- get_legend(
  # create some space to the left of the legend
  CN.T2.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))

## combine plots
isoCNplot<- plot_grid(CvN.T1.mod.plot + theme(legend.position = "none"), 
                    CvN.T2.mod.plot + theme(legend.position = "none"),
                    CN.T1.mod.plot + theme(legend.position = "none"),
                    CN.T2.mod.plot + theme(legend.position = "none"),
                    extract.legend.C.N, rel_widths = c(8,8,8,8,3), ncol=5)

ggsave("figures/isoCNplots.pdf", height=6, width=12, encod="MacRoman")

## combine model diffs
isoCNplot.diff<- plot_grid(CN.diff.T1 + theme(legend.position = "none"), 
                    CN.diff.T2 + theme(legend.position = "none"),
                    rel_widths = c(8,8), ncol=2)

ggsave("figures/isoCN.moddiff.pdf", height=5, width=8, encod="MacRoman")

The “controls” are the tin blanks, plankton, and starting materials

########## ########## ########## 
# run some stats to see how the control material differs from each other

# make a burning treatment
topes.controls$Burn.Trt=as.factor(word(topes.controls$Type, -1, sep="[.]"))

# make a plant name
topes.controls$Plant=as.factor(word(topes.controls$Type, 1, sep="[.]"))

### test some models on controls
# remove blanks and plankton, keeping only plants
topes.controls.plants<-topes.controls[!(topes.controls$Plant=="blank" |
                                          topes.controls$Plant=="plankton" ),]
topes.controls.plants$Plant<-droplevels(topes.controls.plants$Plant)

# keep only non-enriched samples (remove sage)
topes.controls.non.enrich<-topes.controls[!(topes.controls$Plant=="blank" |
                                          topes.controls$Plant=="sage" ),]
topes.controls.non.enrich$Plant<-droplevels(topes.controls.non.enrich$Plant)


######## test plant species differences
mwu(topes.controls.plants, d15N, Plant)
## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = sage (n = 18) | 2 = willow (n = 8):
##   U = 315.000, W = 144.000, p < .001, Z = 4.000
##   effect-size r =  0.784
##    rank-mean(1) = 17.50
##    rank-mean(2) =  4.50
# d15N sage and willow differ (p<0.001)
mwu(topes.controls.plants, C.N, Plant) 
## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = sage (n = 16) | 2 = willow (n = 8):
##   U = 160.000, W = 24.000, p = 0.014, Z = -2.449
##   effect-size r =  0.500
##    rank-mean(1) = 10.00
##    rank-mean(2) = 17.50
# C.N sage and willow differ (p=0.013)

par(mfrow=c(1,2))
boxplot(d15N~Plant, data=topes.controls.plants)
boxplot(C.N~Plant, data=topes.controls.plants)

######## test difference between willow and plankton 
mwu(topes.controls.non.enrich, d15N, Plant)
## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = plankton (n = 7) | 2 = willow (n = 8):
##   U = 28.000, W = 0.000, p = 0.001, Z = -3.240
##   effect-size r =  0.837
##    rank-mean(1) =  4.00
##    rank-mean(2) = 11.50
 # d15N plankton and willow differ (p<0.001)
mwu(topes.controls.non.enrich, C.N, Plant) # C:N plankton and willow differ (p<0.001)
## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = plankton (n = 7) | 2 = willow (n = 8):
##   U = 28.000, W = 0.000, p = 0.001, Z = -3.240
##   effect-size r =  0.837
##    rank-mean(1) =  4.00
##    rank-mean(2) = 11.50
par(mfrow=c(1,2))
boxplot(d15N~Plant, data=topes.controls.non.enrich)
boxplot(C.N~Plant, data=topes.controls.non.enrich)

############# separate plant dfs
#### Sage d15N
topes.controls.sage<-topes.controls[(topes.controls$Plant=="sage"),]
topes.controls.sage$Plant<-droplevels(topes.controls.sage$Plant)
topes.controls.sage$Burn.Trt<-droplevels(topes.controls.sage$Burn.Trt)

# how do different types of sage compare across burn/unburn
# first, no difference between burned or very burned sage
anova(lm(d15N~Burn.Trt, data=topes.controls.sage[!(topes.controls.sage$Burn.Trt=="unburned"),]))
## Analysis of Variance Table
## 
## Response: d15N
##           Df Sum Sq Mean Sq F value Pr(>F)
## Burn.Trt   1   4078  4077.6  1.1836  0.298
## Residuals 12  41339  3444.9
# convert to just 2 levels, no difference here either
topes.controls.sage$Burn.Unb<-ifelse(topes.controls.sage$Burn.Trt=="burned", "burned",
                                     ifelse(topes.controls.sage$Burn.Trt=="very burned", "burned",
                                     "unburned"))


mod.sage<-lm(d15N~Burn.Trt, data=topes.controls.sage) # keep at 3 levels
anova(mod.sage) 
## Analysis of Variance Table
## 
## Response: d15N
##           Df Sum Sq Mean Sq F value Pr(>F)
## Burn.Trt   2   5178  2589.1  0.9125 0.4227
## Residuals 15  42561  2837.4
# no difference in d15N for burned, unburned, very burned sage (p=0.423)

#### Sage C.N
mod.sage.CN<-lm(C.N~Burn.Trt, data=topes.controls.sage)
anova(mod.sage.CN) 
## Analysis of Variance Table
## 
## Response: C.N
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## Burn.Trt   2 3010.5 1505.26   11.32 0.001423 **
## Residuals 13 1728.7  132.98                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(mod.sage.CN, ~Burn.Trt)
multcomp::cld(posthoc, Letters=letters)
##  Burn.Trt   emmean   SE df lower.CL upper.CL .group
##  burned       37.5 4.08 13     28.7     46.3  a    
##  unburned     43.8 5.77 13     31.3     56.2  a    
##  veryburned   70.7 5.77 13     58.2     83.1   b   
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping letter,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
# Sage: difference in unburned, burned, very burned for C:N

par(mfrow=c(1,2))
boxplot(d15N~Burn.Trt, data=topes.controls.sage)
boxplot(C.N~Burn.Trt, data=topes.controls.sage)

###########################
#### Willow d15N
topes.controls.will<-topes.controls[(topes.controls$Plant=="willow"),]
topes.controls.will$Plant<-droplevels(topes.controls.will$Plant)
topes.controls.will$Burn.Trt<-droplevels(topes.controls.will$Burn.Trt)

mod<-lm(C.N~Burn.Trt, data=topes.controls.will)
anova(mod) 
## Analysis of Variance Table
## 
## Response: C.N
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## Burn.Trt   1 14.773 14.7729  5.2794 0.06129 .
## Residuals  6 16.789  2.7982                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# no difference in burned/unburned willow d15N

#### Willow C.N
mod<-lm(C.N~Burn.Trt, data=topes.controls.will)
anova(mod) 
## Analysis of Variance Table
## 
## Response: C.N
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## Burn.Trt   1 14.773 14.7729  5.2794 0.06129 .
## Residuals  6 16.789  2.7982                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Willow: no difference in unburned, burned C:N (p=0.061)

par(mfrow=c(1,2))
boxplot(d15N~Burn.Trt, data=topes.controls.will)
boxplot(C.N~Burn.Trt, data=topes.controls.will)

######### make summary dfs
# summarize by plants
d15N.plant<-aggregate(d15N~Plant, topes.controls, FUN=mean)
d15N.plantSD<-aggregate(d15N~Plant, topes.controls, FUN=length)
d15N.plant[3]<-d15N.plantSD[2]
colnames(d15N.plant)<-c("Plant", "d15N", "SD")

CN.plant<-aggregate(C.N~Plant, topes.controls, FUN=mean)
CN.plantSD<-aggregate(C.N~Plant, topes.controls, FUN=sd)
CN.plant[3]<-CN.plantSD[2]
colnames(CN.plant)<-c("Plant", "C.N", "SD")

# summary df d15N
d15N.cont<-aggregate(d15N~Type, topes.controls, FUN=mean)
d15N.cont.n<-aggregate(d15N~Type, topes.controls, FUN=length)
d15N.cont.SD<-aggregate(d15N~Type, topes.controls, FUN=sd)
d15N.cont[3]<- d15N.cont.SD[2]
d15N.cont[4]<- d15N.cont.n[2]
colnames(d15N.cont)<-c("Type", "d15N", "SD", "n")

# summary df control C:N
CN.cont<-aggregate(C.N~Type, topes.controls, FUN=mean)
CN.cont.n<-aggregate(C.N~Type, topes.controls, FUN=length)
CN.cont.SD<-aggregate(C.N~Type, topes.controls, FUN=sd)
CN.cont[3]<- CN.cont.SD[2]
CN.cont[4]<- CN.cont.n[2]
colnames(CN.cont)<-c("Type", "C.N", "SD", "n")

make boxplots of control sample d15N and C:N

########## control plots
# set levels
topes.controls$Type<-factor(topes.controls$Type, 
                         levels=c("blank", "plankton.stock", 
                                  "willow.unburned", "willow.burned",
                                  "sage.unburned", "sage.burned", 
                                  "sage.veryburned", "sage.stem.burned"))

#### controls d15N boxplot
iso.plot.control.d15N<-ggplot(data=topes.controls,  aes(x=Type, y=d15N, fill=Type)) +
  geom_boxplot() +
  geom_point(pch = 21, position = position_jitterdodge(), alpha=0.6)+
  ylab(expression(paste(delta^{15}, N, " (\u2030, air)"))) +
  scale_fill_manual(values = c("azure2", "cornflowerblue",
                                "darkgoldenrod1", "indianred2",
                                "aquamarine3", "antiquewhite3",
                                "darkolivegreen4", "lightsalmon")) +
  xlab("control types") +  Fig.formatting +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 


###### control C:N
topes.controls.CN<-topes.controls %>% drop_na(C.N) # drop the NAs for C.N, makes plotting problematic

iso.plot.control.CN<-ggplot(data=topes.controls.CN,  aes(x=Type, y=C.N, fill=Type)) +
  geom_boxplot() +
  geom_point(pch = 21, position = position_jitterdodge(), alpha=0.6)+
  ylab("C:N") +
  scale_fill_manual(values = c("cornflowerblue",
                                "darkgoldenrod1", "indianred2",
                                "aquamarine3", "antiquewhite3",
                                "darkolivegreen4", "lightsalmon")) +
  xlab("control types") +  Fig.formatting +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))


####### combine plots
# legend
extract.legend.cont <- get_legend(
  # create some space to the left of the legend
  iso.plot.control.d15N + theme(legend.box.margin = margin(0, 0, 0, 10)))


## combine 
control.iso.alltime<- 
  plot_grid(iso.plot.control.d15N + theme(legend.position = "none"),
            iso.plot.control.CN + theme(legend.position = "none"),
            extract.legend.cont, rel_widths = c(8,8,3), ncol=3)

ggsave("figures/iso.controls.pdf", encod="MacRoman", height=5, width=10)

Use a mixing model to calculate % sage for plankton

# mixing model
head(topes.trt)
topes.trt<-droplevels(topes.trt)


### values for controls
d15N.cont # summary mean d15N by all controls
d15N.plant # summary by plants, # stock plantkon 11, # sage 296, #  willow 13

# summary mean atom percent enrichment
F.cont<-aggregate(at.P..15N  ~ Type, topes.controls, mean)
F.plant<-aggregate(at.P..15N~Plant, topes.controls, FUN=mean)
# sage ~0.475
# willow 0.371

# 2 source mixing model (Post 2002), used d15N values here

# alpha = percent Sage from food web 1
# %Sage = (d15N sample - d15N base 2 [i.e., no-label food])/ (d15N sage food 1 - source 2)
#  d15N values of base 2 = 11 permil for algae/plankton stock/willow
#  d15N value of base 1 = 298 permil for sage

# framed differently from Robinson 2001, TREE
# xtracer = frction of tracer
# Xtracer = (d15N-sample - d15N background) / (d15N-tracer - d15N-background)

topes.trt$percent.sage<-(topes.trt$d15N-12)/(296-12)*100

# unicode text for micrometer = \u03BC, use this in legend

model for % sage

############# all plankton T1
m1.T1.sage <- gam(percent.sage ~ Treatment + Type +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m2.T1.sage <- gam(percent.sage ~ Treatment + Type +
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m3.T1.sage <- gam(percent.sage ~
                     s(plant.mass..g), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")


AIC.sage.T1<-AIC(m1.T1.sage, m2.T1.sage, m3.T1.sage)
## additive model best fit

summary(m1.T1.sage)
anova.gam(m1.T1.sage)
gam.check(m1.T1.sage, rep=1000)
draw(m1.T1.sage)
concrvity(m1.T1.sage)
par(mfrow = c(1, 2))
plot(m1.T1.sage, all.terms = TRUE, page=1)
# model predictions
per.Sage.diff.T1<-plot_difference(
  m1.T1.sage,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
# model for smoothing
msmooth.T1.sage<- gam(percent.sage ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

#plot for the model output on rawdata
per.Sage.T1.mod.plot<-
  plot_smooths(
  model = msmooth.T1.sage,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=percent.sage, color=Treatment, shape=Type)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab("% Sage")+
  xlab("plant material (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(0, 100)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

# overall an effect of burning with both smoothers being significant by treatment
# no effect of type, POM and plankton with similar d15N values


############# all plankton T2
m1.T2.sage <- gam(percent.sage ~ Treatment + Type +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m2.T2.sage <- gam(percent.sage ~ Treatment + Type +
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m3.T2.sage <- gam(percent.sage ~
                     s(plant.mass..g), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")


AIC.sage.T2<-AIC(m1.T2.sage, m2.T2.sage, m3.T2.sage)
# m2.T2.sage preferred
# compare across

summary(m1.T2.sage)
anova.gam(m1.T2.sage)
gam.check(m1.T2.sage, rep=1000)
draw(m1.T2.sage)
concrvity(m1.T2.sage)
par(mfrow = c(1, 2))
plot(m1.T2.sage, all.terms = TRUE, page=1)
# model predictions
per.Sage.diff.T2<-plot_difference(
  m1.T2.sage,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
# model for smoothing
msmooth.T2.sage<- gam(percent.sage ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

#plot for the model output on rawdata
per.Sage.T2.mod.plot<-
  plot_smooths(
  model = msmooth.T2.sage,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=percent.sage, color=Treatment, shape=Type)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab("% Sage")+
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 100)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

# effect of treatment (p<0.001) but not Type (p=0.321)
# smoother significant for both terms

model plots for percent sage mixing model

## AIC table
mod.sag.top<-rep(c( "Treatment + Type + s(plant.mass..g, by=Treatment)", 
              "Treatment + Type + s(plant.mass..g)",
              "s(plant.mass..g)"), times=2)

mod.sag.df<- data.frame(mod.sag.top)
AIC.sag.topes<-bind_rows(AIC.sage.T1, AIC.sage.T2)

AIC.sag.mod<-cbind(mod.sag.df, AIC.sag.topes)
write.csv(AIC.sag.mod, "output/AIC models/AIC.sag.mod.csv")


# legend
extract.legend.mix <- get_legend(
  # create some space to the left of the legend
  per.Sage.T2.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))


## combine 
sage.mix.model<- plot_grid(per.Sage.T1.mod.plot + theme(legend.position = "none"), 
                    per.Sage.T2.mod.plot + theme(legend.position = "none"),
                    extract.legend.mix, rel_widths = c(8,8,3), ncol=3)

ggsave("figures/Isotope.mixmodel.pdf", encod="MacRoman", height=4, width=8)


# and plot difference
sage.mix.model.diff<- plot_grid(per.Sage.diff.T1 + theme(legend.position = "none"), 
                    per.Sage.diff.T2 + theme(legend.position = "none"),
                    extract.legend.mix, rel_widths = c(8,8,3), ncol=3)

ggsave("figures/Isotope.mix.plotdiff.pdf", encod="MacRoman", height=4, width=8)

make d15N plots as well – these follow the % sage, but are informative with the control plots to see the d15N of plankton and the 2 end members.

#### d15N isotope plots for treatments

############# all plankton T1
m1.T1.d15N <- gam(d15N ~ Treatment + Type +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m2.T1.d15N <- gam(d15N ~ Treatment + Type +
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m3.T1.d15N <- gam(d15N ~
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

AIC.d15N.T1<-AIC(m1.T1.d15N, m2.T1.d15N, m3.T1.d15N)
## additive model best fit

summary(m1.T1.d15N)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d15N ~ Treatment + Type + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         96.337      2.610  36.908  < 2e-16 ***
## Treatmentunburned   16.852      3.014   5.591 9.64e-07 ***
## TypePOM              3.954      3.014   1.312    0.196    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value    
## s(plant.mass..g):Treatmentburned   3.560  4.338 136.8  <2e-16 ***
## s(plant.mass..g):Treatmentunburned 3.921  4.762 173.6  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.961   Deviance explained = 96.7%
## -REML = 227.82  Scale est. = 136.26    n = 60
anova.gam(m1.T1.d15N)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d15N ~ Treatment + Type + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df      F  p-value
## Treatment  1 31.261 9.64e-07
## Type       1  1.721    0.196
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   3.560  4.338 136.8  <2e-16
## s(plant.mass..g):Treatmentunburned 3.921  4.762 173.6  <2e-16
gam.check(m1.T1.d15N, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-2.454513e-09,1.329425e-10]
## (score 227.823 & scale 136.2639).
## Hessian positive definite, eigenvalue range [0.5347593,27.64001].
## Model rank =  21 / 21 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 3.56       1    0.47
## s(plant.mass..g):Treatmentunburned 9.00 3.92       1    0.38
draw(m1.T1.d15N)

concrvity(m1.T1.d15N)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 6.67e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.24e-26
## 3 worst    s(plant.mass..g):Treatmentunburned   4.28e-26
## 4 observed para                                 6.67e- 1
## 5 observed s(plant.mass..g):Treatmentburned     7.44e-31
## 6 observed s(plant.mass..g):Treatmentunburned   3.91e-30
## 7 estimate para                                 6.67e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.27e-28
## 9 estimate s(plant.mass..g):Treatmentunburned   1.26e-28
par(mfrow = c(1, 2))
plot(m1.T1.d15N, all.terms = TRUE, page=1)

# model predictions
d15N.diff.T1<-plot_difference(
  m1.T1.d15N,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
# model for smoothing
msmooth.T1.d15N<- gam(d15N ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

#plot for the model output on rawdata
d15N.T1.mod.plot<-
  plot_smooths(
  model = msmooth.T1.d15N,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=d15N, color=Treatment, shape=Type)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste(delta^{15}, N, " (\u2030, air)"))) +
  xlab("plant material (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(0, 250)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

# overall an effect of burning with both smoothers being significant by treatment
# no effect of type, POM and plankton with similar d15N values


############# all plankton T2
m1.T2.d15N <- gam(d15N ~ Treatment + Type +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m2.T2.d15N <- gam(d15N ~ Treatment + Type +
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m3.T2.d15N <- gam(d15N ~
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")


AIC.d15N.T2<-AIC(m1.T2.d15N, m2.T2.d15N, m3.T2.d15N)
## additive model best fit

summary(m1.T2.d15N)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d15N ~ Treatment + Type + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         91.147      2.944  30.960  < 2e-16 ***
## Treatmentunburned   12.290      3.398   3.617 0.000704 ***
## TypePOM             -3.409      3.402  -1.002 0.321237    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F p-value    
## s(plant.mass..g):Treatmentburned   4.669  5.621  79.45  <2e-16 ***
## s(plant.mass..g):Treatmentunburned 3.340  4.082 125.45  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.943   Deviance explained = 95.3%
## -REML = 235.35  Scale est. = 172.72    n = 60
anova.gam(m1.T2.d15N)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d15N ~ Treatment + Type + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df      F  p-value
## Treatment  1 13.082 0.000704
## Type       1  1.004 0.321237
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F p-value
## s(plant.mass..g):Treatmentburned   4.669  5.621  79.45  <2e-16
## s(plant.mass..g):Treatmentunburned 3.340  4.082 125.45  <2e-16
gam.check(m1.T2.d15N, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-7.936647e-08,5.277139e-09]
## (score 235.3485 & scale 172.7225).
## Hessian positive definite, eigenvalue range [1.026571,27.67886].
## Model rank =  21 / 21 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 4.67    0.95    0.32
## s(plant.mass..g):Treatmentunburned 9.00 3.34    0.95    0.28
draw(m1.T2.d15N)

concrvity(m1.T2.d15N)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 0.669   
## 2 worst    s(plant.mass..g):Treatmentburned     0.00464 
## 3 worst    s(plant.mass..g):Treatmentunburned   0.0148  
## 4 observed para                                 0.669   
## 5 observed s(plant.mass..g):Treatmentburned     0.00100 
## 6 observed s(plant.mass..g):Treatmentunburned   0.00307 
## 7 estimate para                                 0.669   
## 8 estimate s(plant.mass..g):Treatmentburned     0.000543
## 9 estimate s(plant.mass..g):Treatmentunburned   0.00169
par(mfrow = c(1, 2))
plot(m1.T2.d15N, all.terms = TRUE, page=1)

# model predictions
d15N.diff.T2<-plot_difference(
  m1.T2.d15N,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
# model for smoothing
msmooth.T2.d15N<- gam(d15N ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

#plot for the model output on rawdata
d15N.T2.mod.plot<-
  plot_smooths(
  model = msmooth.T2.d15N,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=d15N, color=Treatment, shape=Type)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste(delta^{15}, N, " (\u2030, air)"))) +
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 250)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

d15N.model<- plot_grid(d15N.T1.mod.plot + theme(legend.position = "none"), 
                    d15N.T2.mod.plot + theme(legend.position = "none"),
                    extract.legend.mix, rel_widths = c(8,8,3), ncol=3)

ggsave("figures/d15N.model.pdf", encod="MacRoman", height=5, width=10)

################

mod.d15N<-rep(c("Treatment + Type + s(plant.mass..g, by=Treatment)", 
              "Treatment + Type + s(plant.mass..g)",
              "s(plant.mass..g)"), times=2)

mod.d15N.df<- data.frame(mod.d15N)
AIC.d15N<-bind_rows(AIC.d15N.T1, AIC.d15N.T2)

AIC.d15N.mod<-cbind(mod.d15N.df, AIC.d15N)
write.csv(AIC.d15N.mod, "output/AIC models/AIC.d15N.mod.csv")
#################
### 
POM<-topes.trt[(topes.trt$Type=="POM"),]
plankton<-topes.trt[(topes.trt$Type=="plankton"),]

topes <- POM %>% inR_join( plankton, 
                              by=c('Time.point','Treatment', 'plant.mass..g', 'Tank'))


names(topes)[names(topes) == 'percent.sage.x'] <- 'POM.per.sage'
names(topes)[names(topes) == 'percent.sage.y'] <- 'plank.per.sage'

topes<- topes %>% 
  select(Time.point, Treatment, plant.mass..g, Tank, POM.per.sage, plank.per.sage)

topes$TTE.diff<- topes$POM.per.sage - topes$plank.per.sage
topes$TTE.per <- (topes$plank.per.sage / topes$POM.per.sage)*100

topes$TTE.per[topes$TTE.per < 0] <-"" # remove outlier
topes$TTE.per<-as.numeric(topes$TTE.per)

# and again...
topes$TTE.per[topes$TTE.per > 200] <-"" # remove outlier
topes$TTE.per<-as.numeric(topes$TTE.per)



################### TTE as difference
## run some models
m1.T1.TTE <- gam(TTE.diff ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes, method="REML", family="gaussian")

m2.T1.TTE <- gam(TTE.diff ~ Treatment +
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes, method="REML", family="gaussian")

m3.T1.TTE <- gam(TTE.diff ~
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes, method="REML", family="gaussian")


AIC.TTE.T1<-AIC(m1.T1.TTE, m2.T1.TTE, m3.T1.TTE)
## first model best fit

summary(m1.T1.TTE)
anova.gam(m1.T1.TTE)
gam.check(m1.T1.TTE, rep=1000)
draw(m1.T1.TTE)
concrvity(m1.T1.TTE)
par(mfrow = c(1, 2))
plot(m1.T1.TTE, all.terms = TRUE, page=1)

# model predictions
TTE.T1<-plot_difference(
  m1.T1.TTE,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)


###########  
# model for smoothing
msmooth.T1.TTE<- gam(TTE.diff ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes, method="REML", family="gaussian")

#plot for the model output on rawdata
TTE.T1.mod.plot<-
  plot_smooths(
  model = msmooth.T1.TTE,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes[(topes$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=TTE.diff, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab("TTE= POM %sage - Plank %sage")+
  xlab("plant material (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(-10, 25)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))


################### TTE as difference TIME 2
## run some models
m1.T2.TTE <- gam(TTE.diff ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes, method="REML", family="gaussian")

m2.T2.TTE <- gam(TTE.diff ~ Treatment +
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes, method="REML", family="gaussian")

m3.T2.TTE <- gam(TTE.diff ~
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes, method="REML", family="gaussian")


AIC.TTE.T2<-AIC(m1.T2.TTE, m2.T2.TTE, m3.T2.TTE)
## first model best fit

summary(m1.T2.TTE)
anova.gam(m1.T2.TTE)
gam.check(m1.T2.TTE, rep=1000)
draw(m1.T2.TTE)
concrvity(m1.T2.TTE)
par(mfrow = c(1, 2))
plot(m1.T2.TTE, all.terms = TRUE, page=1)

# model predictions
TTE.T2<-plot_difference(
  m1.T2.TTE,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)


###########  
# model for smoothing
msmooth.T2.TTE<- gam(TTE.diff ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes, method="REML", family="gaussian")

#plot for the model output on rawdata
TTE.T2.mod.plot<-
  plot_smooths(
  model = msmooth.T2.TTE,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes[(topes$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=TTE.diff, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab("TTE= POM %sage - Plank %sage")+
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(-10, 25)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))





############ as percent
m1.T1.TTE.per <- gam(TTE.per ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes, method="REML", family="gaussian")

m2.T1.TTE.per <- gam(TTE.per ~ Treatment +
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes, method="REML", family="gaussian")

m3.T1.TTE.per <- gam(TTE.per ~
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes, method="REML", family="gaussian")


AIC.TTE.per.T1<-AIC(m1.T1.TTE.per, m2.T1.TTE.per, m3.T1.TTE.per)
## first model best fit

summary(m1.T1.TTE.per)
anova.gam(m1.T1.TTE.per)
gam.check(m1.T1.TTE.per, rep=1000)
draw(m1.T1.TTE.per)
concrvity(m1.T1.TTE.per)
par(mfrow = c(1, 2))
plot(m1.T1.TTE.per, all.terms = TRUE, page=1)

# model predictions
TTE.per.T1<-plot_difference(
  m1.T1.TTE.per,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)


###########  
# model for smoothing
msmooth.T1.TTE.per<- gam(TTE.per ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes, method="REML", family="gaussian")

#plot for the model output on rawdata
TTE.per.T1.mod.plot<-
  plot_smooths(
  model = msmooth.T1.TTE.per,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes[(topes$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=TTE.per, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab("TTE= (plank - POM) *100")+
  xlab("plant material (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(0, 200)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

### Time 2

m1.T2.TTE.per <- gam(TTE.per ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes, method="REML", family="gaussian")

m2.T2.TTE.per <- gam(TTE.per ~ Treatment +
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes, method="REML", family="gaussian")

m3.T2.TTE.per <- gam(TTE.per ~
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes, method="REML", family="gaussian")


AIC.TTE.per.T2<-AIC(m1.T2.TTE.per, m2.T2.TTE.per, m3.T2.TTE.per)
## first model best fit

summary(m1.T2.TTE.per)
anova.gam(m1.T2.TTE.per)
gam.check(m1.T2.TTE.per, rep=1000)
draw(m1.T2.TTE.per)
concrvity(m1.T2.TTE.per)
par(mfrow = c(1, 2))
plot(m1.T2.TTE.per, all.terms = TRUE, page=1)

# model predictions
TTE.per.T2<-plot_difference(
  m1.T2.TTE.per,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)


###########  
# model for smoothing
msmooth.T2.TTE.per<- gam(TTE.per ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes, method="REML", family="gaussian")

#plot for the model output on rawdata
TTE.per.T2.mod.plot<-
  plot_smooths(
  model = msmooth.T2.TTE.per,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes[(topes$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=TTE.per, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab("TTE= (plank - POM) *100")+
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 200)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

## combine 

TTE.diff.plot<- plot_grid(TTE.T1.mod.plot + theme(legend.position = "none"), 
                    TTE.T2.mod.plot + theme(legend.position = "none"),
                    extract.legend, rel_widths = c(8,8,3), ncol=3)


TTE.per.plot<- plot_grid(TTE.per.T1.mod.plot + theme(legend.position = "none"), 
                    TTE.per.T2.mod.plot + theme(legend.position = "none"),
                    extract.legend, rel_widths = c(8,8,3), ncol=3)

ggsave("figures/TTE.diff.plot.pdf", encod="MacRoman", height=4, width=8)
ggsave("figures/TTE.per.plot.pdf", encod="MacRoman", height=4, width=8)

Plant material

Plant material for the starting material (sage or willow, stems or sticks). This is useful in determining how fire impacted the nutrients in the plant material.

First, run some stats to see what is happening and where differences lie.

plant.nut<-read.csv("data/Pyro_plant material_elemental.csv")

fac<-c("Sample.Name", "type", "plant", "treatment", "rep") # make factors
plant.nut[fac]<-lapply(plant.nut[fac],factor)

sage.nut<-plant.nut[(plant.nut$plant=="sage"),]
will.nut<-plant.nut[(plant.nut$plant=="willow"),]

### Sage ### 
######### %N
plant.N.sag<-lm(N~treatment*type, data=sage.nut)
Anova(plant.N.sag, type=3) # 2 way interaction and main effects
## Anova Table (Type III tests)
## 
## Response: N
##                Sum Sq Df   F value    Pr(>F)    
## (Intercept)    7.0227  1 2171.5213 4.971e-11 ***
## treatment      0.0280  1    8.6632  0.018615 *  
## type           0.2604  1   80.5246 1.894e-05 ***
## treatment:type 0.0705  1   21.8099  0.001602 ** 
## Residuals      0.0259  8                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(plant.N.sag, ~treatment| type)
multcomp::cld(posthoc, Letters=letters)
## type = leaf:
##  treatment emmean     SE df lower.CL upper.CL .group
##  burned     1.530 0.0328  8    1.454     1.61  a    
##  unburned   1.667 0.0328  8    1.591     1.74   b   
## 
## type = stem:
##  treatment emmean     SE df lower.CL upper.CL .group
##  unburned   0.943 0.0328  8    0.868     1.02  a    
##  burned     1.113 0.0328  8    1.038     1.19   b   
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping letter,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
########## %P 
plant.P.sag<-lm(P~treatment*type, data=sage.nut)
Anova(plant.P.sag, type=3) # just type
## Anova Table (Type III tests)
## 
## Response: P
##                  Sum Sq Df  F value    Pr(>F)    
## (Intercept)    0.314928  1 804.0715 2.586e-09 ***
## treatment      0.000771  1   1.9677  0.198292    
## type           0.007561  1  19.3060  0.002306 ** 
## treatment:type 0.000056  1   0.1438  0.714371    
## Residuals      0.003133  8                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
########## %K
plant.K.sag<-lm(K~treatment*type, data=sage.nut)
Anova(plant.K.sag, type=3) # type and treatment
## Anova Table (Type III tests)
## 
## Response: K
##                 Sum Sq Df  F value    Pr(>F)    
## (Intercept)    12.5256  1 372.2328 5.404e-08 ***
## treatment       0.2563  1   7.6157   0.02469 *  
## type            0.1803  1   5.3571   0.04934 *  
## treatment:type  0.0736  1   2.1882   0.17733    
## Residuals       0.2692  8                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(plant.K.sag, ~treatment)
multcomp::cld(posthoc, Letters=letters)
##  treatment emmean     SE df lower.CL upper.CL .group
##  unburned    1.61 0.0749  8     1.44     1.79  a    
##  burned      1.87 0.0749  8     1.70     2.04   b   
## 
## Results are averaged over the levels of: type 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping letter,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
########## %S
plant.S.sag<-lm(S~treatment*type, data=sage.nut)
Anova(plant.S.sag, type=3) # type effect
## Anova Table (Type III tests)
## 
## Response: S
##                 Sum Sq Df  F value    Pr(>F)    
## (Intercept)    0.67783  1 859.3665 1.986e-09 ***
## treatment      0.00010  1   0.1321    0.7257    
## type           0.11152  1 141.3891 2.299e-06 ***
## treatment:type 0.00039  1   0.4885    0.5044    
## Residuals      0.00631  8                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############## Zn ppm
plant.ZN.sag<-lm(ZN~treatment*type, data=sage.nut)
Anova(plant.ZN.sag, type=3) #no effect
## Anova Table (Type III tests)
## 
## Response: ZN
##                 Sum Sq Df F value   Pr(>F)   
## (Intercept)    19959.4  1 15.1723 0.004576 **
## treatment       3073.6  1  2.3364 0.164903   
## type            3280.7  1  2.4938 0.152947   
## treatment:type  1236.3  1  0.9398 0.360731   
## Residuals      10524.1  8                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(plant.ZN.sag, ~treatment| type)
multcomp::cld(posthoc, Letters=letters)
## type = leaf:
##  treatment emmean   SE df lower.CL upper.CL .group
##  unburned    36.3 20.9  8    -12.0     84.6  a    
##  burned      81.6 20.9  8     33.3    129.9  a    
## 
## type = stem:
##  treatment emmean   SE df lower.CL upper.CL .group
##  unburned    30.1 20.9  8    -18.2     78.4  a    
##  burned      34.8 20.9  8    -13.5     83.1  a    
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping letter,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
###############################################
###############################################
### Willow ### 
######### %N
plant.N.wil<-lm(N~treatment*type, data=will.nut)
Anova(plant.N.wil, type=3) #  main effects
## Anova Table (Type III tests)
## 
## Response: N
##                Sum Sq Df  F value    Pr(>F)    
## (Intercept)    8.7381  1 644.0658 3.766e-08 ***
## treatment      0.0817  1   6.0194   0.04388 *  
## type           1.4538  1 107.1530 1.703e-05 ***
## treatment:type 0.0436  1   3.2119   0.11620    
## Residuals      0.0950  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
########## %P 
plant.P.wil<-lm(P~treatment*type, data=will.nut)
Anova(plant.P.wil, type=3) # type and treatment
## Anova Table (Type III tests)
## 
## Response: P
##                  Sum Sq Df  F value    Pr(>F)    
## (Intercept)    0.156865  1 447.9429 1.323e-07 ***
## treatment      0.006403  1  18.2834  0.003674 ** 
## type           0.007239  1  20.6703  0.002647 ** 
## treatment:type 0.000880  1   2.5131  0.156921    
## Residuals      0.002451  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(plant.P.wil, ~treatment)
multcomp::cld(posthoc, Letters=letters)
##  treatment emmean      SE df lower.CL upper.CL .group
##  unburned   0.143 0.00764  7    0.125    0.161  a    
##  burned     0.190 0.00854  7    0.170    0.210   b   
## 
## Results are averaged over the levels of: type 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping letter,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
########## %K
plant.K.wil<-lm(K~treatment*type, data=will.nut)
Anova(plant.K.wil, type=3) # type and treatment
## Anova Table (Type III tests)
## 
## Response: K
##                Sum Sq Df  F value    Pr(>F)    
## (Intercept)    4.6128  1 548.0120  6.59e-08 ***
## treatment      0.0676  1   8.0344 0.0252428 *  
## type           0.2640  1  31.3583 0.0008161 ***
## treatment:type 0.0015  1   0.1725 0.6903484    
## Residuals      0.0589  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(plant.K.wil, ~treatment)
multcomp::cld(posthoc, Letters=letters)
##  treatment emmean     SE df lower.CL upper.CL .group
##  unburned   0.817 0.0375  7    0.728    0.905  a    
##  burned     1.006 0.0419  7    0.906    1.105   b   
## 
## Results are averaged over the levels of: type 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping letter,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
########## %S
plant.S.wil<-lm(S~treatment*type, data=will.nut)
Anova(plant.S.wil, type=3) # main and interactions
## Anova Table (Type III tests)
## 
## Response: S
##                  Sum Sq Df  F value    Pr(>F)    
## (Intercept)    0.197633  1 1748.117 1.168e-09 ***
## treatment      0.002521  1   22.303  0.002151 ** 
## type           0.042375  1  374.819 2.446e-07 ***
## treatment:type 0.001355  1   11.985  0.010519 *  
## Residuals      0.000791  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(plant.S.wil, ~treatment| type)
multcomp::cld(posthoc, Letters=letters)
## type = leaf:
##  treatment emmean      SE df lower.CL upper.CL .group
##  unburned  0.2157 0.00614  7   0.2012   0.2302  a    
##  burned    0.2567 0.00614  7   0.2422   0.2712   b   
## 
## type = stem:
##  treatment emmean      SE df lower.CL upper.CL .group
##  burned    0.0688 0.00752  7   0.0510   0.0865  a    
##  unburned  0.0728 0.00614  7   0.0583   0.0873  a    
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping letter,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
############## Zn ppm
plant.ZN.wil<-lm(ZN~treatment*type, data=will.nut)
Anova(plant.ZN.wil, type=3) #no effect
## Anova Table (Type III tests)
## 
## Response: ZN
##                Sum Sq Df  F value    Pr(>F)    
## (Intercept)    124440  1 196.4564 2.229e-06 ***
## treatment       22363  1  35.3043 0.0005748 ***
## type            21956  1  34.6631 0.0006071 ***
## treatment:type   5340  1   8.4306 0.0228697 *  
## Residuals        4434  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(plant.ZN.wil, ~treatment| type)
multcomp::cld(posthoc, Letters=letters)
## type = leaf:
##  treatment emmean   SE df lower.CL upper.CL .group
##  unburned    81.6 14.5  7    47.21    115.9  a    
##  burned     203.7 14.5  7   169.31    238.0   b   
## 
## type = stem:
##  treatment emmean   SE df lower.CL upper.CL .group
##  unburned    35.8 14.5  7     1.44     70.2  a    
##  burned      68.4 17.8  7    26.32    110.5  a    
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping letter,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

Make a figure of the elemental analysis data.

plant.nut$int.fac<-factor(interaction(plant.nut$plant, plant.nut$type),
                         levels=c("sage.leaf", "sage.stem", "willow.leaf", "willow.stem"))
####### figures
N.plot<- ggplot(plant.nut, aes(x=int.fac, y=N, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75))+
  xlab("Plant:Tissue")+
  ylab("%N")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

P.plot<- ggplot(plant.nut, aes(x=int.fac, y=P, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75))+
  xlab("Plant:Tissue")+
  ylab("%P")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

K.plot<- ggplot(plant.nut, aes(x=int.fac, y=K, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75))+
  xlab("Plant:Tissue")+
  ylab("%K")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))


S.plot<- ggplot(plant.nut, aes(x=int.fac, y=S, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75))+
  xlab("Plant:Tissue")+
  ylab("%S")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Zn.plot<- ggplot(plant.nut, aes(x=int.fac, y=ZN, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75)) +
  xlab("Plant:Tissue")+
  ylab("Zn ppm")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
  
extract.legend.nut <- get_legend(
  # create some space to the left of the legend
  Zn.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))

nutrients<-plot_grid(
  N.plot+ theme(legend.position = "none"),
  P.plot+ theme(legend.position = "none"),
  K.plot+ theme(legend.position = "none"), 
  S.plot+ theme(legend.position = "none"), 
  Zn.plot+ theme(legend.position = "none"),
  extract.legend.nut, 
  rel_widths = c(8,8,8,8,8,3), ncol=6)

nutrients

ggsave("figures/leaf.nutrients.pdf", height=4, width=12)

H2O Phosphorus

The phosphorous in water was only run at Time-2. Make a plot here and run the model.

phos<-read.csv("data/Pyro_water.phosph.csv")

fac2<-c("Time.Point", "Treatment", "Tank") # make factors
phos[fac2]<-lapply(phos[fac2],factor)

phos<-na.omit(phos)

############# phosphorous T2: significant smoothers, not treatment overall
m1.T2.phos<- gam(TP.umol..l ~ Treatment +
            s(plant.mass..g, by=Treatment), data=phos, method="REML", family="gaussian")

m2.T2.phos<- gam(TP.umol..l ~ Treatment +
            s(plant.mass..g), data=phos, method="REML", family="gaussian")

m3.T2.phos<- gam(TP.umol..l ~
            s(plant.mass..g), data=phos, method="REML", family="gaussian")

AIC.T2.phos<-AIC(m1.T2.phos,m2.T2.phos, m3.T2.phos)
# factor smooth best

summary(m1.T2.phos)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TP.umol..l ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        10.3649     0.5518  18.785  6.6e-12 ***
## Treatmentunburned  -0.8831     0.7659  -1.153    0.267    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value    
## s(plant.mass..g):Treatmentburned   2.924  3.525 154.7  <2e-16 ***
## s(plant.mass..g):Treatmentunburned 2.930  3.371 124.6  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.979   Deviance explained = 98.5%
## -REML = 46.397  Scale est. = 2.9168    n = 23
anova.gam(m1.T2.phos)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TP.umol..l ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 1.329   0.267
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.924  3.525 154.7  <2e-16
## s(plant.mass..g):Treatmentunburned 2.930  3.371 124.6  <2e-16
gam.check(m1.T2.phos, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-8.634731e-07,9.305738e-08]
## (score 46.39739 & scale 2.916807).
## Hessian positive definite, eigenvalue range [0.3315034,9.7044].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value  
## s(plant.mass..g):Treatmentburned   9.00 2.92    0.75    0.08 .
## s(plant.mass..g):Treatmentunburned 9.00 2.93    0.75    0.07 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
draw(m1.T2.phos)

concrvity(m1.T2.phos)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                   0.662 
## 2 worst    s(plant.mass..g):Treatmentburned       0.293 
## 3 worst    s(plant.mass..g):Treatmentunburned     1.00  
## 4 observed para                                   0.662 
## 5 observed s(plant.mass..g):Treatmentburned       0.0265
## 6 observed s(plant.mass..g):Treatmentunburned     0.0382
## 7 estimate para                                   0.662 
## 8 estimate s(plant.mass..g):Treatmentburned       0.0410
## 9 estimate s(plant.mass..g):Treatmentunburned     0.0314
par(mfrow = c(1, 2))
plot(m1.T2.phos, all.terms = TRUE, page=1)

# model predictions
phos.diff.T2<-plot_difference(
  m1.T2.phos,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
phos.T2.mod.plot<-
  plot_smooths(
  model = m1.T2.phos,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=phos, 
             aes(x=plant.mass..g, y=TP.umol..l, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("TP", ~(mu*mol/L), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))


phos.plots<-plot_grid(phos.T2.mod.plot, phos.diff.T2, ncol=2)
ggsave("figures/phos.plots.pdf", height=5, width=8)


## AIC table
mod.phos<-rep(c("Treatment + s(plant.mass..g, by=Treatment)", 
              "Treatment + s(plant.mass..g)", 
              "s(plant.mass..g)"), times=1)

mod.phos.df<- data.frame(mod.phos)

AIC.phos.mod<-cbind(mod.phos.df, AIC.T2.phos)
write.csv(AIC.phos.mod, "output/AIC models/AIC.phos.mod.csv")

Greenhouse Gas Data

Import data and clean up

GHG<-read.csv("data/GH.gases/Pyro_ghg_headspace.csv")
GHG<-na.omit(GHG) # remove NAs

# set structure
make.fac<-c("Time.point", "Treatment", "Tank", "Gas")
GHG[make.fac] <- lapply(GHG[make.fac], factor) # make all these factors
GHG$plant.mass..g<-as.numeric(GHG$plant.mass..g)

# rename
GHG<- GHG %>% rename("GHG.nM" = "GHGwater..nmol.GHG..L.H2O")

# reorder columns and export
GHG<- GHG %>% select(Time.point, Treatment, plant.mass..g, Tank, Gas, GHG.nM)

GHG<-GHG[!(GHG$GHG.nM < 0),]

Make plots and run models for CO2

CO2<-subset(GHG, Gas=="CO2")

#convert GHG units from nM to uM
CO2$GHG.uM<- CO2$GHG.nM/1000

############# GHG plots and modesl
#### CO2
#######-- T0

m1.T0.CO2<- gam(GHG.uM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T0", data = CO2,
            method="REML", family="gaussian")

m2.T0.CO2<- gam(GHG.uM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T0", data = CO2,
            method="REML", family="gaussian")

m3.T0.CO2<- gam(GHG.uM ~
            s(plant.mass..g), subset = Time.point=="T0", data = CO2,
            method="REML", family="gaussian")

CO2.T0.AIC<-AIC(m1.T0.CO2,m2.T0.CO2, m3.T0.CO2)
#factor smooth best

summary(m1.T0.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        34.4608     0.7814  44.099  < 2e-16 ***
## Treatmentunburned  -4.2776     1.1051  -3.871 0.000776 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F p-value   
## s(plant.mass..g):Treatmentburned   1.000   1.00 10.566 0.00353 **
## s(plant.mass..g):Treatmentunburned 4.022   4.88  3.383 0.02381 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.573   Deviance explained = 66.1%
## -REML = 74.819  Scale est. = 9.1598    n = 30
anova.gam(m1.T0.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F  p-value
## Treatment  1 14.98 0.000776
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F p-value
## s(plant.mass..g):Treatmentburned   1.000  1.000 10.566 0.00353
## s(plant.mass..g):Treatmentunburned 4.022  4.880  3.383 0.02381
gam.check(m1.T0.CO2, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-1.834502e-05,7.70351e-07]
## (score 74.81935 & scale 9.159824).
## Hessian positive definite, eigenvalue range [1.834448e-05,13.18306].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 1.00    1.03    0.48
## s(plant.mass..g):Treatmentunburned 9.00 4.02    1.03    0.47
draw(m1.T0.CO2)

concrvity(m1.T0.CO2)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5.0 e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     1.98e-26
## 3 worst    s(plant.mass..g):Treatmentunburned   2.09e-26
## 4 observed para                                 5.0 e- 1
## 5 observed s(plant.mass..g):Treatmentburned     9.97e-30
## 6 observed s(plant.mass..g):Treatmentunburned   3.55e-29
## 7 estimate para                                 5.00e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     7.95e-29
## 9 estimate s(plant.mass..g):Treatmentunburned   7.49e-29
par(mfrow = c(1, 2))
plot(m1.T0.CO2, all.terms = TRUE, page=1)

# model predictions
CO2.diff.T0<-plot_difference(
  m1.T0.CO2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CO2.T0.mod.plot<-
  plot_smooths(
  model = m1.T0.CO2,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CO2[(CO2$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=GHG.uM, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CO"[2], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-0") +
  coord_cartesian(ylim=c(0, 400)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CO2.T0.mod.plot

#### CO2 
#######-- T1

m1.T1.CO2<- gam(GHG.uM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T1", data = CO2,
            method="REML", family="gaussian")

m2.T1.CO2<- gam(GHG.uM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T1", data = CO2,
            method="REML", family="gaussian")

m3.T1.CO2<- gam(GHG.uM ~
            s(plant.mass..g), subset = Time.point=="T1", data = CO2,
            method="REML", family="gaussian")


CO2.T1.AIC<-AIC(m1.T1.CO2,m2.T1.CO2, m3.T1.CO2)
#factor smooth best

summary(m1.T1.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        137.112      4.830  28.387  < 2e-16 ***
## Treatmentunburned   20.946      6.831   3.066  0.00547 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value    
## s(plant.mass..g):Treatmentburned   2.047  2.527 157.7  <2e-16 ***
## s(plant.mass..g):Treatmentunburned 2.966  3.633 144.2  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.97   Deviance explained = 97.6%
## -REML = 121.27  Scale est. = 349.95    n = 30
anova.gam(m1.T1.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 9.403 0.00547
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.047  2.527 157.7  <2e-16
## s(plant.mass..g):Treatmentunburned 2.966  3.633 144.2  <2e-16
gam.check(m1.T1.CO2, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-1.86084e-07,1.75551e-07]
## (score 121.2725 & scale 349.9528).
## Hessian positive definite, eigenvalue range [0.1792501,13.09977].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 2.05    1.02    0.52
## s(plant.mass..g):Treatmentunburned 9.00 2.97    1.02    0.55
draw(m1.T1.CO2)

concrvity(m1.T1.CO2)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5.0 e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     1.98e-26
## 3 worst    s(plant.mass..g):Treatmentunburned   2.09e-26
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     1.01e-29
## 6 observed s(plant.mass..g):Treatmentunburned   2.32e-30
## 7 estimate para                                 5.00e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     7.95e-29
## 9 estimate s(plant.mass..g):Treatmentunburned   7.49e-29
par(mfrow = c(1, 2))
plot(m1.T1.CO2, all.terms = TRUE, page=1)

# model predictions
CO2.diff.T1<-plot_difference(
  m1.T1.CO2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CO2.T1.mod.plot<-
  plot_smooths(
  model = m1.T1.CO2,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CO2[(CO2$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=GHG.uM, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CO"[2], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(0, 400)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CO2.T1.mod.plot

#### CO2
#######-- T2

m1.T2.CO2<- gam(GHG.uM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T2", data = CO2,
            method="REML", family="gaussian")

m2.T2.CO2<- gam(GHG.uM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T2", data = CO2,
            method="REML", family="gaussian")

m3.T2.CO2<- gam(GHG.uM ~
            s(plant.mass..g), subset = Time.point=="T2", data = CO2,
            method="REML", family="gaussian")

CO2.T2.AIC<-AIC(m1.T2.CO2, m2.T2.CO2, m3.T2.CO2)
# smooth by factor best

summary(m1.T2.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         69.394      4.628  14.995 4.65e-13 ***
## Treatmentunburned    4.275      6.545   0.653     0.52    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value    
## s(plant.mass..g):Treatmentburned   2.499  3.072 52.20  <2e-16 ***
## s(plant.mass..g):Treatmentunburned 3.422  4.176 22.47  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.897   Deviance explained = 92.2%
## -REML = 121.41  Scale est. = 321.23    n = 30
anova.gam(m1.T2.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.427    0.52
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.499  3.072 52.20  <2e-16
## s(plant.mass..g):Treatmentunburned 3.422  4.176 22.47  <2e-16
gam.check(m1.T2.CO2, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 3 iterations.
## Gradient range [-9.205725e-07,1.593428e-06]
## (score 121.4081 & scale 321.2311).
## Hessian positive definite, eigenvalue range [0.3528807,13.16313].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 2.50    0.93    0.26
## s(plant.mass..g):Treatmentunburned 9.00 3.42    0.93    0.33
draw(m1.T2.CO2)

concrvity(m1.T2.CO2)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5.0 e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     1.98e-26
## 3 worst    s(plant.mass..g):Treatmentunburned   2.09e-26
## 4 observed para                                 5.0 e- 1
## 5 observed s(plant.mass..g):Treatmentburned     8.86e-30
## 6 observed s(plant.mass..g):Treatmentunburned   3.18e-30
## 7 estimate para                                 5.00e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     7.95e-29
## 9 estimate s(plant.mass..g):Treatmentunburned   7.49e-29
par(mfrow = c(1, 2))
plot(m1.T2.CO2, all.terms = TRUE, page=1)

# model predictions
CO2.diff.T2<-plot_difference(
  m1.T2.CO2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CO2.T2.mod.plot<-
  plot_smooths(
  model = m1.T2.CO2,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CO2[(CO2$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=GHG.uM, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CO"[2], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 400)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CO2.T2.mod.plot

#### CO2 
#######-- T3

m1.T3.CO2<- gam(GHG.uM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T3", data = CO2,
            method="REML", family="gaussian")

m2.T3.CO2<- gam(GHG.uM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T3", data = CO2,
            method="REML", family="gaussian")

m3.T3.CO2<- gam(GHG.uM ~
            s(plant.mass..g), subset = Time.point=="T3", data = CO2,
            method="REML", family="gaussian")

CO2.T3.AIC<-AIC(m1.T3.CO2, m2.T3.CO2, m3.T3.CO2)
# smooth by factor best

summary(m1.T3.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         36.020      5.742   6.273 2.02e-06 ***
## Treatmentunburned   12.209      7.980   1.530     0.14    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value    
## s(plant.mass..g):Treatmentburned   2.744  3.366 11.86 4.97e-05 ***
## s(plant.mass..g):Treatmentunburned 1.000  1.000 28.09 2.26e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.705   Deviance explained = 75.5%
## -REML = 119.18  Scale est. = 460.52    n = 29
anova.gam(m1.T3.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 2.341    0.14
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value
## s(plant.mass..g):Treatmentburned   2.744  3.366 11.86 4.97e-05
## s(plant.mass..g):Treatmentunburned 1.000  1.000 28.09 2.26e-05
gam.check(m1.T3.CO2, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 9 iterations.
## Gradient range [-3.53958e-05,8.858783e-05]
## (score 119.1836 & scale 460.5193).
## Hessian positive definite, eigenvalue range [3.540667e-05,12.56498].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 2.74    0.96    0.32
## s(plant.mass..g):Treatmentunburned 9.00 1.00    0.96    0.40
draw(m1.T3.CO2)

concrvity(m1.T3.CO2)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 0.536   
## 2 worst    s(plant.mass..g):Treatmentburned     0.0380  
## 3 worst    s(plant.mass..g):Treatmentunburned   0.0113  
## 4 observed para                                 0.536   
## 5 observed s(plant.mass..g):Treatmentburned     0.000735
## 6 observed s(plant.mass..g):Treatmentunburned   0.000138
## 7 estimate para                                 0.536   
## 8 estimate s(plant.mass..g):Treatmentburned     0.000272
## 9 estimate s(plant.mass..g):Treatmentunburned   0.000250
par(mfrow = c(1, 2))
plot(m1.T3.CO2, all.terms = TRUE, page=1)

# model predictions
CO2.diff.T3<-plot_difference(
  m1.T3.CO2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CO2.T3.mod.plot<-
  plot_smooths(
  model = m1.T3.CO2,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CO2[(CO2$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=GHG.uM, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CO"[2], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-3") +
  coord_cartesian(ylim=c(0, 400)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CO2.T3.mod.plot

# models and raw data
CO2.model.plots<-plot_grid(
  CO2.T0.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-0"),
  CO2.T1.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-10"),
  CO2.T2.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-31"),
  CO2.T3.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-59"),
  extract.legend, 
  rel_widths = c(8,8,8,8,3), ncol=5)

### model differences
CO2.plot.diff<-plot_grid(
  CO2.diff.T0+ ggtitle("CO2-T0"),
  CO2.diff.T1+ ggtitle("Day-10"),
  CO2.diff.T2+ ggtitle("Day-31"),
  CO2.diff.T3+ ggtitle("Day-59"),
  rel_widths = c(8,8,8,8), ncol=4)

## AIC table
mod.ghg<-rep(c("Treatment +s(plant.mass..g, by=Treatment)", 
              "Treatment + s(plant.mass..g)", 
              "s(plant.mass..g)"), times=4)

mod.ghg.df<- data.frame(mod.ghg)

AIC.CO2<-bind_rows(CO2.T0.AIC, CO2.T1.AIC, CO2.T2.AIC, CO2.T3.AIC)
AIC.CO2.mod<-cbind(mod.ghg.df, AIC.CO2)
write.csv(AIC.CO2.mod, "output/AIC models/AIC.CO2.mod.csv")

Make plots and run models for methane.

CH4<-subset(GHG, Gas=="CH4")


m1.T0.CH4<- gam(GHG.nM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T0", data = CH4,
            method="REML", family="gaussian")

m2.T0.CH4<- gam(GHG.nM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T0", data = CH4,
            method="REML", family="gaussian")

m3.T0.CH4<- gam(GHG.nM ~
            s(plant.mass..g), subset = Time.point=="T0", data = CH4,
            method="REML", family="gaussian")

CH4.T0.AIC<-AIC(m1.T0.CH4, m2.T0.CH4, m3.T0.CH4)
#factor single smooth best

summary(m2.T0.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.nM ~ Treatment + s(plant.mass..g)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         4.3461     0.6747   6.441 1.16e-06 ***
## Treatmentunburned   1.3989     0.9799   1.427    0.166    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf Ref.df     F p-value
## s(plant.mass..g)   1      1 0.718   0.405
## 
## R-sq.(adj) =  0.0155   Deviance explained = 9.13%
## -REML = 60.338  Scale est. = 6.2802    n = 27
anova.gam(m2.T0.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.nM ~ Treatment + s(plant.mass..g)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 2.038   0.166
## 
## Approximate significance of smooth terms:
##                  edf Ref.df     F p-value
## s(plant.mass..g)   1      1 0.718   0.405
gam.check(m2.T0.CH4, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-1.631119e-05,8.871995e-05]
## (score 60.33814 & scale 6.280154).
## Hessian positive definite, eigenvalue range [1.632563e-05,11.99991].
## Model rank =  11 / 11 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                  k' edf k-index p-value
## s(plant.mass..g)  9   1    1.13    0.68
draw(m2.T0.CH4)

concrvity(m2.T0.CH4)
## # A tibble: 6 × 3
##   type     term             concurvity
##   <chr>    <chr>                 <dbl>
## 1 worst    para                 0.499 
## 2 worst    s(plant.mass..g)     0.0659
## 3 observed para                 0.499 
## 4 observed s(plant.mass..g)     0.0298
## 5 estimate para                 0.499 
## 6 estimate s(plant.mass..g)     0.0232
par(mfrow = c(1, 2))
plot(m2.T0.CH4, all.terms = TRUE, page=1)

# model predictions
CH4.diff.T0<-plot_difference(
  m2.T0.CH4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CH4.T0.mod.plot<-
  plot_smooths(
  model = m2.T0.CH4,
  series = plant.mass..g,
  comparison=Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CH4[(CH4$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=GHG.nM, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  ylab(expression(paste("CH"[4], ~(nM), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-0") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CH4.T0.mod.plot

#### CH4 
#######-- T1

m1.T1.CH4<- gam(GHG.nM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T1", data = CH4,
            method="REML", family="gaussian")

m2.T1.CH4<- gam(GHG.nM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T1", data = CH4,
            method="REML", family="gaussian")

m3.T1.CH4<- gam(GHG.nM ~
            s(plant.mass..g), subset = Time.point=="T1", data = CH4,
            method="REML", family="gaussian")

CH4.T1.AIC<-AIC(m1.T1.CH4, m2.T1.CH4, m3.T1.CH4)
#factor by smooth best

summary(m1.T1.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.nM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        11.6731     0.7459  15.649 1.07e-13 ***
## Treatmentunburned   0.5438     1.0549   0.516    0.611    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df    F p-value   
## s(plant.mass..g):Treatmentburned   1.813  2.244 0.89 0.42743   
## s(plant.mass..g):Treatmentunburned 3.346  4.086 6.53 0.00119 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.476   Deviance explained = 58.7%
## -REML = 73.049  Scale est. = 8.3465    n = 30
anova.gam(m1.T1.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.nM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.266   0.611
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df    F p-value
## s(plant.mass..g):Treatmentburned   1.813  2.244 0.89 0.42743
## s(plant.mass..g):Treatmentunburned 3.346  4.086 6.53 0.00119
gam.check(m1.T1.CH4, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-2.015212e-07,2.357495e-08]
## (score 73.04861 & scale 8.346463).
## Hessian positive definite, eigenvalue range [0.2435469,13.121].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 1.81    1.04    0.56
## s(plant.mass..g):Treatmentunburned 9.00 3.35    1.04    0.56
draw(m1.T1.CH4)

concrvity(m1.T1.CH4)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5.0 e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     1.98e-26
## 3 worst    s(plant.mass..g):Treatmentunburned   2.09e-26
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     1.11e-30
## 6 observed s(plant.mass..g):Treatmentunburned   4.58e-30
## 7 estimate para                                 5.00e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     7.95e-29
## 9 estimate s(plant.mass..g):Treatmentunburned   7.49e-29
par(mfrow = c(1, 2))
plot(m1.T1.CH4, all.terms = TRUE, page=1)

# model predictions
CH4.diff.T1<-plot_difference(
  m1.T1.CH4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CH4.T1.mod.plot<-
  plot_smooths(
  model = m1.T1.CH4,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CH4[(CH4$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=GHG.nM, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CH"[4], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CH4.T1.mod.plot

#### CH4
#######-- T2

m1.T2.CH4<- gam(GHG.nM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T2", data = CH4,
            method="REML", family="gaussian")

m2.T2.CH4<- gam(GHG.nM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T2", data = CH4,
            method="REML", family="gaussian")

m3.T2.CH4<- gam(GHG.nM ~
            s(plant.mass..g), subset = Time.point=="T2", data = CH4,
            method="REML", family="gaussian")

CH4.T2.AIC<-AIC(m1.T2.CH4, m2.T2.CH4, m3.T2.CH4)
# global best

summary(m3.T2.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.nM ~ s(plant.mass..g)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.628      1.066   9.031 8.69e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf Ref.df    F p-value
## s(plant.mass..g)   1  1.001 0.19   0.667
## 
## R-sq.(adj) =  -0.0287   Deviance explained = 0.675%
## -REML = 92.541  Scale est. = 34.099    n = 30
anova.gam(m3.T2.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.nM ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F p-value
## s(plant.mass..g) 1.000  1.001 0.19   0.667
gam.check(m3.T2.CH4, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-2.147014e-05,0.000126711]
## (score 92.54116 & scale 34.09853).
## Hessian positive definite, eigenvalue range [2.148928e-05,13.99987].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                  k' edf k-index p-value
## s(plant.mass..g)  9   1    0.84    0.16
draw(m3.T2.CH4)

concrvity(m3.T2.CH4)
## # A tibble: 6 × 3
##   type     term             concurvity
##   <chr>    <chr>                 <dbl>
## 1 worst    para               2.08e-26
## 2 worst    s(plant.mass..g)   2.09e-26
## 3 observed para               2.08e-26
## 4 observed s(plant.mass..g)   2.62e-32
## 5 estimate para               2.08e-26
## 6 estimate s(plant.mass..g)   7.15e-29
par(mfrow = c(1, 2))
plot(m3.T2.CH4, all.terms = TRUE, page=1)

# model predictions
CH4.diff.T2<-plot_difference(
  m1.T2.CH4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CH4.T2.mod.plot<-
  plot_smooths(
  model = m3.T2.CH4,
  series = plant.mass..g,
)  + theme(legend.position = "none") +
  geom_point(data=CH4[(CH4$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=GHG.nM, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CH"[4], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CH4.T2.mod.plot

#### CH4 
#######-- T3

m1.T3.CH4<- gam(GHG.nM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T3", data = CH4,
            method="REML", family="gaussian")

m2.T3.CH4<- gam(GHG.nM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T3", data = CH4,
            method="REML", family="gaussian")

m3.T3.CH4<- gam(GHG.nM ~
            s(plant.mass..g), subset = Time.point=="T3", data = CH4,
            method="REML", family="gaussian")

CH4.T3.AIC<-AIC(m1.T3.CH4, m2.T3.CH4, m3.T3.CH4)
# global with treatment best

summary(m2.T3.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.nM ~ Treatment + s(plant.mass..g)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         22.395      2.255   9.930 4.33e-10 ***
## Treatmentunburned   -6.089      3.189  -1.909    0.068 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df     F p-value
## s(plant.mass..g) 3.381  4.127 2.038   0.113
## 
## R-sq.(adj) =  0.266   Deviance explained = 37.7%
## -REML = 103.88  Scale est. = 76.293    n = 30
anova.gam(m2.T3.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GHG.nM ~ Treatment + s(plant.mass..g)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 3.645   0.068
## 
## Approximate significance of smooth terms:
##                    edf Ref.df     F p-value
## s(plant.mass..g) 3.381  4.127 2.038   0.113
gam.check(m2.T3.CH4, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-1.041682e-09,9.185275e-11]
## (score 103.8784 & scale 76.2932).
## Hessian positive definite, eigenvalue range [0.4424777,13.60852].
## Model rank =  11 / 11 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value
## s(plant.mass..g) 9.00 3.38    1.15     0.8
draw(m2.T3.CH4)

concrvity(m2.T3.CH4)
## # A tibble: 6 × 3
##   type     term             concurvity
##   <chr>    <chr>                 <dbl>
## 1 worst    para               5   e- 1
## 2 worst    s(plant.mass..g)   2.09e-26
## 3 observed para               5   e- 1
## 4 observed s(plant.mass..g)   5.75e-31
## 5 estimate para               5.0 e- 1
## 6 estimate s(plant.mass..g)   7.15e-29
par(mfrow = c(1, 2))
plot(m2.T3.CH4, all.terms = TRUE, page=1)

# model predictions
CH4.diff.T3<-plot_difference(
  m1.T3.CH4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CH4.T3.mod.plot<-
  plot_smooths(
  model = m2.T3.CH4,
  series = plant.mass..g,
  comparison=Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CH4[(CH4$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=GHG.nM, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  ylab(expression(paste("CH"[4], ~(nM), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-3") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))


CH4.model.plots<-plot_grid(
  CH4.T0.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-0"),
  CH4.T1.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-10"),
  CH4.T2.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-31"),
  CH4.T3.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-59"), 
  extract.legend, 
  rel_widths = c(8,8,8,8,3), ncol=5)

### model differences
CH4.plot.diff<-plot_grid(
  CH4.diff.T0+ ggtitle("CH4-Day-0"),
  CH4.diff.T1+ ggtitle("Day-10"),
  CH4.diff.T2+ ggtitle("Day-31"),
  CH4.diff.T3+ ggtitle("Day-59"),
  rel_widths = c(8,8,8,8), ncol=4)


##### AIC table
AIC.CH4<-bind_rows(CH4.T0.AIC, CH4.T1.AIC, CH4.T2.AIC, CH4.T3.AIC)

AIC.CH4.mod<-cbind(mod.ghg.df, AIC.CH4)
write.csv(AIC.CH4.mod, "output/AIC models/AIC.CH4.mod.csv")

Compile the greenhouse gas plots

GHG.plots<-plot_grid(CO2.model.plots, CH4.model.plots, ncol=1)
ggsave("figures/GHG.molar.mod.plots.pdf", height=7, width=12)


GHG.plot.diff<-plot_grid(CO2.plot.diff, CH4.plot.diff, ncol=1)
ggsave("figures/GHG.molar.mod.diff.pdf", height=7, width=12)